Skip to content

cylindra.components

CylSpline

Bases: Spline

A spline object with cylindrical structure.

Source code in cylindra/components/spline/_cyl_spline.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
class CylSpline(Spline):
    """A spline object with cylindrical structure."""

    @property
    def radius(self) -> nm | None:
        """Average radius of the cylinder."""
        return self.props.get_glob(H.radius, None)

    @radius.setter
    def radius(self, value: nm | None):
        if value is None:
            if H.radius in self.props.glob.columns:
                self.props.drop_glob(H.radius)
            return None
        value = float(value)
        if value <= 0:
            raise ValueError("Radius must be positive.")
        elif value != value:
            raise ValueError("Radius got NaN value.")
        self.props.update_glob([pl.Series(H.radius, [value], dtype=pl.Float32)])
        return None

    def radius_range(self, rc: nm | None = None) -> tuple[nm, nm]:
        """Return the range of the radius used for the cylindric coordinate."""
        if rc is None:
            if self.radius is None:
                raise ValueError("Radius is not set.")
            rc = self.radius
        cfg = self.config
        return (max(rc - cfg.thickness_inner, 0.0), rc + cfg.thickness_outer)

    @property
    def orientation(self) -> Ori:
        """Orientation of the spline."""
        return Ori(str(self.props.get_glob(H.orientation, "none")))

    @orientation.setter
    def orientation(self, value: Ori | str | None):
        if value is None:
            value = Ori.none
        else:
            value = Ori(value)
        self.props.update_glob([pl.Series(H.orientation, [str(value)])])
        return None

    def invert(self) -> CylSpline:
        """
        Invert the direction of spline. Also invert orientation if exists.

        Returns
        -------
        CylSpline
            Inverted object
        """
        # NOTE: invert() calls clip() internally.
        # We don't have to invert the orientation here.
        new = super().invert()
        new.props.update_loc(self.props.loc[::-1], self.props.window_size)
        return new

    def clip(self, start: float, stop: float) -> CylSpline:
        """
        Clip spline and generate a new one.

        This method does not convert spline bases. ``_lims`` is updated instead.
        For instance, if you want to clip spline at 20% to 80% position, call
        ``spl.clip(0.2, 0.8)``. If ``stop < start``, the orientation of spline
        will be inverted, thus the ``orientation`` attribute will also be inverted.

        Parameters
        ----------
        start : float
            New starting position.
        stop : float
            New stopping position.

        Returns
        -------
        CylSpline
            Clipped spline.
        """
        clipped = super().clip(start, stop)

        clipped.props.glob = self.props.glob.clone()
        clipped.props._binsize_glob = self.props._binsize_glob.copy()
        if start > stop:
            clipped.orientation = Ori.invert(self.orientation)
        else:
            clipped.orientation = self.orientation
        return clipped

    def restore(self) -> CylSpline:
        """
        Restore the original, not-clipped spline.

        Returns
        -------
        Spline
            Copy of the original spline.
        """
        original = super().restore()
        start, stop = self._lims
        if start > stop:
            original.orientation = Ori.invert(self.orientation)
        else:
            original.orientation = self.orientation
        return original

    def nrise(self, **kwargs) -> int:
        """Raw start number (minus for microtubule)"""
        if H.start in kwargs:
            start = kwargs[H.start]
        elif H.start in self.props.glob.columns:
            start = self.props.get_glob(H.start)
        else:
            start = self.cylinder_params(**kwargs).start
        return start * self.config.rise_sign

    def cylinder_params(self, **kwargs) -> CylinderParameters:
        """Get the cylinder parameters of the spline."""
        radius = _get_globalprops(self, kwargs, H.radius)
        if radius is None:
            raise ValueError("Radius is not known.")
        radius += (self.config.thickness_outer - self.config.thickness_inner) / 2
        return CylinderParameters.solve(
            spacing=_get_globalprops(self, kwargs, H.spacing),
            pitch=_get_globalprops(self, kwargs, H.pitch),
            twist=_get_globalprops(self, kwargs, H.twist),
            skew=_get_globalprops(self, kwargs, H.skew),
            rise_angle=_get_globalprops(self, kwargs, H.rise),
            radius=radius,
            npf=_get_globalprops(self, kwargs, H.npf),
            start=_get_globalprops(self, kwargs, H.start),
            allow_duplicate=True,
            rise_sign=self.config.rise_sign,
        )

    def cylinder_model(
        self,
        offsets: tuple[float, float] = (0.0, 0.0),
        **kwargs,
    ) -> CylinderModel:
        """
        Return the cylinder model of the spline.

        Parameters
        ----------
        offsets : tuple of float, optional
            Offset of the model. See :meth:`map_monomers` for details.

        Returns
        -------
        CylinderModel
            The cylinder model.
        """
        length = self.length()
        cp = self.cylinder_params(**kwargs)
        ly = cp.spacing_proj
        la = cp.lat_spacing_proj
        factor = ly / la
        ny = roundint(length / ly) + 1  # number of monomers in y-direction

        if offsets is None:
            offsets = (0.0, 0.0)

        return CylinderModel(
            shape=(ny, cp.npf),
            tilts=(
                cp.tan_skew * factor,
                cp.tan_rise_raw / factor,
            ),
            intervals=(ly, la / cp.perimeter * 2 * np.pi),
            radius=cp.radius,
            offsets=offsets,
        )

    def update_props(
        self,
        *,
        npf: int | None = None,
        start: int | None = None,
        orientation: Ori | str | None = None,
    ) -> CylSpline:
        """Update the npf or orientation parameters in place."""
        loc = list[pl.Expr]()
        glob = list[pl.Series]()
        if npf is not None:
            loc.append(pl.repeat(npf, pl.len()).cast(pl.UInt8).alias(H.npf))
            glob.append(pl.Series([npf]).cast(pl.UInt8).alias(H.npf))
        if start is not None:
            loc.append(pl.repeat(start, pl.len()).cast(pl.Int8).alias(H.start))
            glob.append(pl.Series([start]).cast(pl.Int8).alias(H.start))
        if orientation is not None:
            glob.append(
                pl.Series([str(orientation)]).cast(pl.String).alias(H.orientation)
            )

        self.props.loc = with_columns(self.props.loc, loc)
        self.props.glob = with_columns(self.props.glob, glob)

        return self

    def update_glob_by_cylinder_params(self, cparams: CylinderParameters) -> CylSpline:
        """Update the global properties using a CylinderParameters object."""
        props = {
            H.rise: cparams.rise_angle,
            H.rise_length: cparams.rise_length,
            H.pitch: cparams.pitch,
            H.spacing: cparams.spacing,
            H.skew: cparams.skew,
            H.twist: cparams.twist,
            H.npf: cparams.npf,
            H.start: cparams.start,
            H.radius: cparams.radius,
        }
        self.props.update_glob(props)
        return self

    def _need_rotation(self, orientation: Ori | str | None) -> bool:
        if orientation is not None:
            orientation = Ori(orientation)
            if orientation is Ori.none or self.orientation is Ori.none:
                raise ValueError(
                    "Either molecules' orientation or the input orientation should "
                    "not be none."
                )
            if orientation is not self.orientation:
                return True
        return False
orientation property writable

Orientation of the spline.

radius property writable

Average radius of the cylinder.

clip(start, stop)

Clip spline and generate a new one.

This method does not convert spline bases. _lims is updated instead. For instance, if you want to clip spline at 20% to 80% position, call spl.clip(0.2, 0.8). If stop < start, the orientation of spline will be inverted, thus the orientation attribute will also be inverted.

Parameters:

Name Type Description Default
start float

New starting position.

 required
stop float

New stopping position.

 required

Returns:

Type Description
CylSpline

Clipped spline.
Source code in cylindra/components/spline/_cyl_spline.py 
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def clip(self, start: float, stop: float) -> CylSpline:
    """
    Clip spline and generate a new one.

    This method does not convert spline bases. ``_lims`` is updated instead.
    For instance, if you want to clip spline at 20% to 80% position, call
    ``spl.clip(0.2, 0.8)``. If ``stop < start``, the orientation of spline
    will be inverted, thus the ``orientation`` attribute will also be inverted.

    Parameters
    ----------
    start : float
        New starting position.
    stop : float
        New stopping position.

    Returns
    -------
    CylSpline
        Clipped spline.
    """
    clipped = super().clip(start, stop)

    clipped.props.glob = self.props.glob.clone()
    clipped.props._binsize_glob = self.props._binsize_glob.copy()
    if start > stop:
        clipped.orientation = Ori.invert(self.orientation)
    else:
        clipped.orientation = self.orientation
    return clipped
cylinder_model(offsets=(0.0, 0.0), **kwargs)

Return the cylinder model of the spline.

Parameters:

Name Type Description Default
offsets tuple of float

Offset of the model. See :meth:map_monomers for details.

 (0.0, 0.0)

Returns:

Type Description
CylinderModel

The cylinder model.
Source code in cylindra/components/spline/_cyl_spline.py 
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
def cylinder_model(
    self,
    offsets: tuple[float, float] = (0.0, 0.0),
    **kwargs,
) -> CylinderModel:
    """
    Return the cylinder model of the spline.

    Parameters
    ----------
    offsets : tuple of float, optional
        Offset of the model. See :meth:`map_monomers` for details.

    Returns
    -------
    CylinderModel
        The cylinder model.
    """
    length = self.length()
    cp = self.cylinder_params(**kwargs)
    ly = cp.spacing_proj
    la = cp.lat_spacing_proj
    factor = ly / la
    ny = roundint(length / ly) + 1  # number of monomers in y-direction

    if offsets is None:
        offsets = (0.0, 0.0)

    return CylinderModel(
        shape=(ny, cp.npf),
        tilts=(
            cp.tan_skew * factor,
            cp.tan_rise_raw / factor,
        ),
        intervals=(ly, la / cp.perimeter * 2 * np.pi),
        radius=cp.radius,
        offsets=offsets,
    )
cylinder_params(**kwargs)

Get the cylinder parameters of the spline.

Source code in cylindra/components/spline/_cyl_spline.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def cylinder_params(self, **kwargs) -> CylinderParameters:
    """Get the cylinder parameters of the spline."""
    radius = _get_globalprops(self, kwargs, H.radius)
    if radius is None:
        raise ValueError("Radius is not known.")
    radius += (self.config.thickness_outer - self.config.thickness_inner) / 2
    return CylinderParameters.solve(
        spacing=_get_globalprops(self, kwargs, H.spacing),
        pitch=_get_globalprops(self, kwargs, H.pitch),
        twist=_get_globalprops(self, kwargs, H.twist),
        skew=_get_globalprops(self, kwargs, H.skew),
        rise_angle=_get_globalprops(self, kwargs, H.rise),
        radius=radius,
        npf=_get_globalprops(self, kwargs, H.npf),
        start=_get_globalprops(self, kwargs, H.start),
        allow_duplicate=True,
        rise_sign=self.config.rise_sign,
    )
invert()

Invert the direction of spline. Also invert orientation if exists.

Returns:

Type Description
CylSpline

Inverted object
Source code in cylindra/components/spline/_cyl_spline.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def invert(self) -> CylSpline:
    """
    Invert the direction of spline. Also invert orientation if exists.

    Returns
    -------
    CylSpline
        Inverted object
    """
    # NOTE: invert() calls clip() internally.
    # We don't have to invert the orientation here.
    new = super().invert()
    new.props.update_loc(self.props.loc[::-1], self.props.window_size)
    return new
nrise(**kwargs)

Raw start number (minus for microtubule)

Source code in cylindra/components/spline/_cyl_spline.py 
124
125
126
127
128
129
130
131
132
def nrise(self, **kwargs) -> int:
    """Raw start number (minus for microtubule)"""
    if H.start in kwargs:
        start = kwargs[H.start]
    elif H.start in self.props.glob.columns:
        start = self.props.get_glob(H.start)
    else:
        start = self.cylinder_params(**kwargs).start
    return start * self.config.rise_sign
radius_range(rc=None)

Return the range of the radius used for the cylindric coordinate.

Source code in cylindra/components/spline/_cyl_spline.py 
38
39
40
41
42
43
44
45
def radius_range(self, rc: nm | None = None) -> tuple[nm, nm]:
    """Return the range of the radius used for the cylindric coordinate."""
    if rc is None:
        if self.radius is None:
            raise ValueError("Radius is not set.")
        rc = self.radius
    cfg = self.config
    return (max(rc - cfg.thickness_inner, 0.0), rc + cfg.thickness_outer)
restore()

Restore the original, not-clipped spline.

Returns:

Type Description
Spline

Copy of the original spline.
Source code in cylindra/components/spline/_cyl_spline.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def restore(self) -> CylSpline:
    """
    Restore the original, not-clipped spline.

    Returns
    -------
    Spline
        Copy of the original spline.
    """
    original = super().restore()
    start, stop = self._lims
    if start > stop:
        original.orientation = Ori.invert(self.orientation)
    else:
        original.orientation = self.orientation
    return original
update_glob_by_cylinder_params(cparams)

Update the global properties using a CylinderParameters object.

Source code in cylindra/components/spline/_cyl_spline.py 
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def update_glob_by_cylinder_params(self, cparams: CylinderParameters) -> CylSpline:
    """Update the global properties using a CylinderParameters object."""
    props = {
        H.rise: cparams.rise_angle,
        H.rise_length: cparams.rise_length,
        H.pitch: cparams.pitch,
        H.spacing: cparams.spacing,
        H.skew: cparams.skew,
        H.twist: cparams.twist,
        H.npf: cparams.npf,
        H.start: cparams.start,
        H.radius: cparams.radius,
    }
    self.props.update_glob(props)
    return self
update_props(*, npf=None, start=None, orientation=None)

Update the npf or orientation parameters in place.

Source code in cylindra/components/spline/_cyl_spline.py 
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def update_props(
    self,
    *,
    npf: int | None = None,
    start: int | None = None,
    orientation: Ori | str | None = None,
) -> CylSpline:
    """Update the npf or orientation parameters in place."""
    loc = list[pl.Expr]()
    glob = list[pl.Series]()
    if npf is not None:
        loc.append(pl.repeat(npf, pl.len()).cast(pl.UInt8).alias(H.npf))
        glob.append(pl.Series([npf]).cast(pl.UInt8).alias(H.npf))
    if start is not None:
        loc.append(pl.repeat(start, pl.len()).cast(pl.Int8).alias(H.start))
        glob.append(pl.Series([start]).cast(pl.Int8).alias(H.start))
    if orientation is not None:
        glob.append(
            pl.Series([str(orientation)]).cast(pl.String).alias(H.orientation)
        )

    self.props.loc = with_columns(self.props.loc, loc)
    self.props.glob = with_columns(self.props.glob, glob)

    return self

Spline

Bases: BaseComponent

3D spline curve model with coordinate system.

Anchor points can be set via anchor property. A spline object is semi-immutable. Different spline curves are always of different objects, but the anchors and properties can be dynamically changed.

References
  • Scipy document https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html
Source code in cylindra/components/spline/_spline_base.py 
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
class Spline(BaseComponent):
    """
    3D spline curve model with coordinate system.

    Anchor points can be set via `anchor` property. A spline object is semi-immutable.
    Different spline curves are always of different objects, but the anchors and
    properties can be dynamically changed.

    References
    ----------
    - Scipy document
      https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html
    """

    def __init__(
        self,
        order: int = 3,
        *,
        lims: tuple[float, float] = (0.0, 1.0),
        extrapolate: ExtrapolationMode | str = ExtrapolationMode.linear,
        config: dict[str, Any] | SplineConfig = {},
    ):
        self._tck: TCKType = (None, None, order)
        self._u: NDArray[np.float32] | None = None
        self._anchors = None
        self._extrapolate = ExtrapolationMode(extrapolate)

        self._lims = lims
        self._props = SplineProps()
        if isinstance(config, SplineConfig):
            self._config = config
        else:
            self._config = SplineConfig.construct(**config)

    @property
    def props(self) -> SplineProps:
        """Return the spline properties"""
        return self._props

    @property
    def config(self) -> SplineConfig:
        """Return the spline configuration"""
        return self._config

    def has_props(self) -> bool:
        """True if there are any properties."""
        return len(self.props.loc) > 0 or len(self.props.glob) > 0

    def copy(self, copy_props: bool = True, copy_config: bool = True) -> Self:
        """
        Copy Spline object.

        Parameters
        ----------
        copy_props : bool, default True
            Also copy local/global properties if true.

        Returns
        -------
        Spline
            Copied object.
        """
        new = self.__class__(
            order=self.order, lims=self._lims, extrapolate=self.extrapolate
        )
        new._tck = self._tck
        new._u = self._u
        new._anchors = self._anchors

        if copy_props:
            new._props = self.props.copy()
        if copy_config:
            new._config = self.config.copy()
        return new

    __copy__ = copy

    def with_extrapolation(self, extrapolate: ExtrapolationMode | str) -> Self:
        """Return a copy of the spline with a new extrapolation mode."""
        new = self.copy()
        new._extrapolate = ExtrapolationMode(extrapolate)
        return new

    def with_config(
        self,
        config: dict[str, Any] | SplineConfig,
        copy_props: bool = False,
    ) -> Self:
        """Return a copy of the spline with a new config."""
        new = self.copy(copy_props=copy_props, copy_config=False)
        if not isinstance(config, SplineConfig):
            config = SplineConfig.construct(**config)
        new._config = config
        return new

    @property
    def knots(self) -> np.ndarray:
        """Spline knots."""
        return self._tck[0]

    @property
    def coeff(self) -> list[np.ndarray]:
        """Spline coefficient."""
        return self._tck[1]

    @property
    def order(self) -> int:
        """Spline order."""
        return self._tck[2]

    @property
    def extrapolate(self) -> ExtrapolationMode:
        """Extrapolation mode of the spline."""
        return self._extrapolate

    @property
    def params(self) -> np.ndarray:
        """Spline parameters."""
        return self._u

    @classmethod
    def line(cls, start: ArrayLike, end: ArrayLike) -> Self:
        """
        Create a line spline.

        Parameters
        ----------
        start : array_like
            Start point of the line.
        end : array_like
            End point of the line.

        Returns
        -------
        Spline
            Line spline.
        """
        spl = cls()
        coords = np.stack([start, end], axis=0)
        return spl.fit(coords, err_max=0.0)

    def translate(self, shift: tuple[nm, nm, nm]):
        """Translate the spline by given shift vectors."""
        new = self.copy()
        c = [x + s for x, s in zip(self.coeff, shift, strict=True)]
        new._tck = (self.knots, c, self.order)
        return new

    @property
    def has_anchors(self) -> bool:
        """True if there are any anchors."""
        return self._anchors is not None

    @property
    def anchors(self) -> NDArray[np.float32]:
        """Local anchors along spline."""
        if self._anchors is None:
            raise ValueError("Anchor has not been set yet.")
        return self._anchors

    @anchors.setter
    def anchors(self, positions: float | Sequence[float]) -> None:
        _anc = np.atleast_1d(np.asarray(positions, dtype=np.float32))
        if _anc.ndim != 1:
            raise TypeError("Could not convert positions into 1D array.")
        elif _anc.min() < 0 or _anc.max() > 1:
            msg = (
                "Anchor positions should be set between 0 and 1. Otherwise spline "
                "curve does not fit well."
            )
            warnings.warn(msg, UserWarning, stacklevel=2)
        _old = self._anchors
        if (
            _old is None
            or _anc.size != _old.size
            or not np.allclose(_anc, _old, rtol=1e-4, atol=1e-4)
        ):
            self.props.clear_loc()
        self._anchors = _anc
        return None

    @anchors.deleter
    def anchors(self) -> None:
        self._anchors = None
        self.props.clear_loc()
        return None

    def is_inverted(self) -> bool:
        """Return true if spline is inverted."""
        return self._lims[0] > self._lims[1]

    @property
    def lims(self) -> tuple[float, float]:
        """Return spline limit positions."""
        return self._lims

    def _set_params(self, tck, u) -> Self:
        self._tck = tck
        self._u = u
        return self

    def prep_anchor_positions(
        self,
        interval: nm | None = None,
        n: int | None = None,
        max_interval: nm | None = None,
    ) -> NDArray[np.float32]:
        length = self.length()
        if interval is not None:
            stop, n_segs = interval_divmod(length, interval)
            n = n_segs + 1
        elif n is not None:
            stop = length
        elif max_interval is not None:
            n = max(ceilint(length / max_interval), self.order) + 1
            stop = length
        else:
            raise ValueError("Either 'interval' or 'n' must be specified.")
        y = np.linspace(0, stop, n)
        return self.y_to_position(y)

    def make_anchors(
        self,
        interval: nm | None = None,
        n: int | None = None,
        max_interval: nm | None = None,
    ) -> Self:
        """
        Make anchor points at constant intervals. Either interval, number of anchor or
        the maximum interval between anchors can be specified.

        Parameters
        ----------
        interval : nm, optional
            Interval between anchor points.
        n : int, optional
            Number of anchor points, including both ends.
        max_interval: nm, optional
            Spline will be split by as little anchors as possible but interval between
            anchors will not be larger than this. The number of anchors are also
            guaranteed to be larger than spline order.
        """
        self.anchors = self.prep_anchor_positions(interval, n, max_interval)
        return self

    def __repr__(self) -> str:
        """Use start/end points to describe a spline."""
        start, end = self.map(self._lims)
        start = "({:.1f}, {:.1f}, {:.1f})".format(*start)
        end = "({:.1f}, {:.1f}, {:.1f})".format(*end)
        return f"Spline[{start}:{end}]"

    def close_to(self: Self, other: Self) -> bool:
        """True if two objects draws the same curve."""
        if not isinstance(other, self.__class__):
            return False
        t0, c0, k0 = self._tck
        t1, c1, k1 = other._tck
        return (
            np.allclose(t0, t1)
            and all(np.allclose(x, y) for x, y in zip(c0, c1, strict=True))
            and k0 == k1
            and np.allclose(self._u, other._u)
            and np.allclose(self._lims, other._lims)
        )

    def clip(self, start: float, stop: float) -> Self:
        """
        Clip spline and generate a new one.

        This method does not convert spline bases. `_lims` is updated instead.

        Parameters
        ----------
        start : float
            New starting position.
        stop : float
            New stopping position.

        Returns
        -------
        Spline
            Clipped spline.
        """
        u0 = _linear_conversion(start, *self._lims)
        u1 = _linear_conversion(stop, *self._lims)
        return self.__class__(
            order=self.order,
            lims=(u0, u1),
            extrapolate=self.extrapolate,
            config=self.config,
        )._set_params(self._tck, self._u)

    def restore(self) -> Self:
        """
        Restore the original, not-clipped spline.

        Returns
        -------
        Spline
            Copy of the original spline.
        """
        return self.__class__(
            order=self.order,
            lims=(0, 1),
            extrapolate=self.extrapolate,
            config=self.config,
        )._set_params(self._tck, self._u)

    def resample(self, max_interval: nm = 1.0, err_max: nm = 0.1) -> Self:
        """
        Resample a new spline along the original spline.

        Parameters
        ----------
        max_interval : nm, default 1.0
            Maximum interval between resampling points.
        err_max : float, default 0.1
            Spline fitting maximum error.

        Returns
        -------
        Spline
            Resampled spline object.
        """
        l = self.length()
        points = self.map(np.linspace(0, 1, ceilint(l / max_interval)))
        return self.fit(points, err_max=err_max)

    def fit(
        self,
        coords: ArrayLike,
        *,
        err_max: nm = 1.0,
    ) -> Self:
        """
        Fit spline model to coordinates.

        This method uses `scipy.interpolate.splprep` to fit given coordinates to a
        spline.

        Parameters
        ----------
        coords : np.ndarray
            Coordinates. Must be (N, 3).
        err_max : float, default 1.0
            Error allowed for fitting. Several upper limit of residual values will be
            used and the fit that results in error lower than this value and minimize
            the maximum curvature will be chosen.

        Returns
        -------
        Spline
            New spline fit to given coordinates.
        """
        crds = np.asarray(coords)
        npoints = crds.shape[0]
        if npoints < 2:
            raise ValueError("Number of input coordinates must be > 1.")
        if npoints <= self.order:
            k = npoints - 1
        else:
            k = self.order

        if self.is_inverted():
            crds = crds[::-1]

        if err_max > 4.0:
            raise ValueError("std_max must be smaller than 4.0.")
        if err_max < 1e-3:
            std_list = [err_max]
        else:
            ntrial = max(int(err_max / 0.02), 2)
            std_list = np.linspace(0, err_max, ntrial)[1:]

        fit_results = list[SplineFitResult]()
        new = self.__class__(order=k, extrapolate=self.extrapolate, config=self.config)
        with warnings.catch_warnings():
            # fitting may fail for some std
            warnings.simplefilter("ignore", RuntimeWarning)
            for std in std_list:
                _tck, _u = splprep(crds.T, k=k, s=std**2 * npoints)
                new._set_params(_tck, _u)
                _crds_at_u = new.map(_u)
                res: NDArray[np.float32] = np.sqrt(
                    np.sum((_crds_at_u - crds) ** 2, axis=1)
                )
                _knots = _tck[0][new.order : -new.order]
                nedge = _knots.size - 1
                assert nedge > 0
                nanc = nedge * 20 + 1
                anc = np.interp(
                    np.linspace(0, 1, nanc), np.linspace(0, 1, nedge + 1), _knots
                )
                max_curvature = new.curvature(anc).max()
                success = res.max() <= err_max
                fit_results.append(
                    SplineFitResult((_tck, _u), max_curvature, res, success)
                )

        fit_results_filt = list(filter(lambda x: x.success, fit_results))
        if len(fit_results_filt) == 0:
            fit_results_filt = fit_results

        reult_opt = min(fit_results_filt, key=lambda x: x.curvature)
        return new._set_params(*reult_opt.params)

    def shift(
        self,
        positions: Sequence[float] | None = None,
        shifts: NDArray[np.floating] | None = None,
        *,
        err_max: nm = 1.0,
    ) -> Self:
        """
        Fit spline model using a list of shifts in XZ-plane.

        Parameters
        ----------
        positions : sequence of float, optional
            Positions. Between 0 and 1. If not given, anchors are used instead.
        shifts : np.ndarray
            Shift from center in nm. Must be (N, 2).
        err_max : float, default 1.0
            Error allowed for fitting. See `Spline.fit`.

        Returns
        -------
        Spline
            Spline shifted by fitting to given coordinates.
        """
        if shifts is None:
            raise ValueError("Shifts must be given.")
        coords = self.map(positions)
        rot = self.get_rotator(positions)
        # insert 0 in y coordinates.
        shifts = np.stack([shifts[:, 0], np.zeros(len(rot)), shifts[:, 1]], axis=1)
        coords += rot.apply(shifts)
        return self.fit(coords, err_max=err_max)

    def distances(
        self,
        positions: Sequence[float] | None = None,
        nknots: int = 512,
    ) -> NDArray[np.float32]:
        """
        Get the distances from u=0.

        Parameters
        ----------
        positions : sequence of float, optional
            Positions. Between 0 and 1. If not given, anchors are used instead.

        Returns
        -------
        np.ndarray
            Distances for each `u`.
        """
        if positions is None:
            _u = self.anchors
        else:
            _u = np.asarray(positions, dtype=np.float32)
            if _u.ndim != 1:
                raise ValueError("Positions must be 1D array.")
        u = np.linspace(0, 1, nknots)
        u_tr = _linear_conversion(u, *self._lims)
        dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
        dist = np.concatenate([[0], np.sqrt(dx**2 + dy**2 + dz**2)]).cumsum()
        tck = splrep(u, dist, k=1)
        out = splev(_u, tck)
        return out

    def map(
        self,
        positions: float | NDArray[np.number] | None = None,
        der: int = 0,
    ) -> NDArray[np.float32]:
        """
        Calculate coordinates (or n-th derivative) at points on the spline.

        Parameters
        ----------
        positions : np.ndarray or float, optional
            Positions. Between 0 and 1. If not given, anchors are used instead.
        der : int, default 0
            `der`-th derivative will be calculated.

        Returns
        -------
        np.ndarray
            Positions or vectors in (3,) or (N, 3) shape.
        """
        _assert_fitted(self)
        if positions is None:
            positions = self.anchors
        u0, u1 = self._lims
        if np.isscalar(positions):
            if self.order < der:
                return np.zeros(3, dtype=np.float32)
            u_tr = _linear_conversion(float(positions), u0, u1)
            if 0 <= u_tr <= 1 or self.extrapolate is ExtrapolationMode.default:
                coord = splev([u_tr], self._tck, der=der)
            elif self.extrapolate is ExtrapolationMode.linear:
                if der == 0:
                    if u_tr < 0:
                        der0 = splev([0], self._tck, der=0)
                        der1 = splev([0], self._tck, der=1)
                        dr = u_tr
                    else:
                        der0 = splev([1], self._tck, der=0)
                        der1 = splev([1], self._tck, der=1)
                        dr = u_tr - 1
                    coord = [a0 + a1 * dr for a0, a1 in zip(der0, der1, strict=True)]
                elif der == 1:
                    if u_tr < 0:
                        coord = splev([0], self._tck, der=1)
                    else:
                        coord = splev([1], self._tck, der=1)
                else:
                    coord = [[0], [0], [0]]
            else:
                raise ValueError(f"Invalid extrapolation mode: {self.extrapolate!r}.")
            out = np.concatenate(coord).astype(np.float32)

        else:
            u_tr = _linear_conversion(np.asarray(positions, dtype=np.float32), u0, u1)
            if self.order < der:
                return np.zeros((u_tr.size, 3), dtype=np.float32)
            if self.extrapolate is ExtrapolationMode.default:
                out = np.stack(splev(u_tr, self._tck, der=der), axis=1).astype(
                    np.float32
                )
            elif self.extrapolate is ExtrapolationMode.linear:
                sl_small = u_tr < 0
                sl_large = u_tr > 1
                n_small = np.count_nonzero(sl_small)
                n_large = np.count_nonzero(sl_large)
                out = np.stack(splev(u_tr, self._tck, der=der), axis=1).astype(
                    np.float32
                )
                if der == 0:
                    if n_small > 0:
                        der0 = np.array(splev(0, self._tck, der=0), dtype=np.float32)
                        der1 = np.array(splev(0, self._tck, der=1), dtype=np.float32)
                        dr = u_tr[sl_small]
                        coords_new = der0 + der1 * dr[:, np.newaxis]
                        out[sl_small] = coords_new

                    if n_large > 0:
                        der0 = splev(1, self._tck, der=0)
                        der1 = splev(1, self._tck, der=1)
                        dr = u_tr[sl_large] - 1
                        coords_new = der0 + der1 * dr[:, np.newaxis]
                        out[sl_large] = coords_new

                elif der == 1:
                    if n_small > 0:
                        out[sl_small] = splev(0, self._tck, der=1)
                    if n_large > 0:
                        out[sl_large] = splev(1, self._tck, der=1)
                else:
                    if n_small > 0:
                        out[sl_small] = 0
                    if n_large > 0:
                        out[sl_large] = 0
            else:
                raise ValueError(f"Invalid extrapolation mode: {self.extrapolate!r}.")

        if u0 > u1 and der % 2 == 1:
            out = -out
        return out

    __call__ = map  # scipy-style alias

    def partition(self, n: int, der: int = 0) -> NDArray[np.float32]:
        """Return the n-partitioning coordinates of the spline."""
        u = np.linspace(0, 1, n)
        return self.map(u, der)

    def length(self, start: float = 0, stop: float = 1, nknots: int = 512) -> nm:
        """
        Approximate the length of B-spline between [start, stop] by partitioning
        the spline with 'nknots' knots. nknots=256 is large enough for most cases.
        """
        _assert_fitted(self)
        u = np.linspace(start, stop, nknots)
        u_tr = _linear_conversion(u, *self._lims)
        dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
        return np.sum(np.sqrt(dx**2 + dy**2 + dz**2))

    def invert(self) -> Self:
        """
        Invert the direction of spline.

        Returns
        -------
        Spline
            Inverted object
        """
        anchors = self._anchors
        inverted = self.clip(1.0, 0.0)
        if anchors is not None:
            inverted.anchors = 1 - anchors[::-1]
        return inverted

    def split(
        self,
        at: nm,
        from_start: bool = True,
        trim: nm = 0.0,
        allow_discard: bool = False,
    ) -> list[Self]:
        spl_len = self.length()
        if not from_start:
            at = spl_len - at
        at_rel_0 = (at - trim) / spl_len
        at_rel_1 = (at + trim) / spl_len
        out = list["Self"]()
        if at_rel_0 >= 0:
            out.append(self.clip(0.0, at_rel_0))
        elif not allow_discard:
            raise ValueError(
                "Split position must be over `trim` if allow_discard is False, but "
                f"tried to split at {at:.1f} nm."
            )
        if at_rel_1 <= 1.0:
            out.append(self.clip(at_rel_1, 1.0))
        elif not allow_discard:
            raise ValueError(
                "Split position must be under `length - trim` if allow_discard is "
                f"False, but tried to split at {at:.1f} nm (length = {spl_len:.1f})."
            )
        return out

    def curvature(
        self,
        positions: Sequence[float] | None = None,
    ) -> NDArray[np.float32]:
        """
        Calculate curvature of spline curve.

        Parameters
        ----------
        positions : sequence of float, optional
            Positions. Between 0 and 1. If not given, anchors are used instead.

        Returns
        -------
        np.ndarray
            Array of curvature.

        References
        ----------
        - https://en.wikipedia.org/wiki/Curvature#Space_curves
        """

        if positions is None:
            positions = self.anchors

        dz, dy, dx = self.map(positions, der=1).T
        ddz, ddy, ddx = self.map(positions, der=2).T
        a = (
            (ddz * dy - ddy * dz) ** 2
            + (ddx * dz - ddz * dx) ** 2
            + (ddy * dx - ddx * dy) ** 2
        )
        return np.sqrt(a) / (dx**2 + dy**2 + dz**2) ** 1.5

    def curvature_radii(self, positions: Sequence[float] = None) -> NDArray[np.float32]:
        """Inverse of curvature."""
        return 1.0 / self.curvature(positions)

    def to_dict(self) -> SplineInfo:
        """Convert spline info into a dict."""
        t, c, k = self._tck
        u = self._u
        return {
            "t": t.tolist(),
            "c": {"z": c[0].tolist(), "y": c[1].tolist(), "x": c[2].tolist()},
            "k": k,
            "u": u.tolist(),
            "lims": self._lims,
            "localprops_window_size": dict(self.props.window_size),
            "binsize_loc": dict(self.props.binsize_loc),
            "binsize_glob": dict(self.props.binsize_glob),
            "extrapolate": self._extrapolate.name,
            "config": self.config.asdict(),
        }

    @classmethod
    def from_dict(cls: type[Self], d: SplineInfo) -> Self:
        """
        Construct a spline model from a dictionary.

        Parameters
        ----------
        d: dict
            Dictionary with keys "t", "c", "k", "u" and "lims".

        Returns
        -------
        Spline
            Spline object constructed from the dictionary.
        """
        self = cls(
            order=d.get("k", 3),
            lims=d.get("lims", (0, 1)),
            extrapolate=d.get("extrapolate", "linear"),
        )
        t = np.asarray(d["t"])
        c = [np.asarray(d["c"][k]) for k in "zyx"]
        k = roundint(d["k"])
        self._tck = (t, c, k)
        self._u = np.asarray(d["u"])
        self.props._window_size = d.get("localprops_window_size", {})
        self.props._binsize_loc = d.get("binsize_loc", {})
        self.props._binsize_glob = d.get("binsize_glob", {})
        if cfg := d.get("config", None):
            self._config = SplineConfig.from_dict(cfg)
        return self

    def get_rotator(
        self,
        positions: Sequence[float] | None = None,
        inverse: bool = False,
    ) -> Rotation:
        """
        Calculate list of Affine transformation matrix along spline, which correspond to
        the orientation of spline curve.

        Parameters
        ----------
        positions : array-like, (N,)
            Positions. Between 0 and 1.
        inverse : bool, default False
            If True, rotation matrix will be inversed.

        Returns
        -------
        Rotation
            Rotation object at each anchor.
        """
        if positions is None:
            positions = self.anchors
        ds = self.map(positions, der=1)
        out = axes_to_rotator(None, -ds)

        if inverse:
            out = out.inv()

        return out

    def local_cartesian(
        self,
        shape: tuple[nm, nm],
        depth: nm,
        u: float | Sequence[float] = None,
        scale: nm = 1.0,
    ) -> NDArray[np.float32]:
        """
        Generate local Cartesian coordinate systems.

        The generated array can be used for `ndi.map_coordinates`. The result coordinate
        systems are flat, i.e., not distorted by the curvature of spline.

        Parameters
        ----------
        shape : (float, float)
            Vertical and horizontal length of Cartesian coordinates. Corresponds to zx
            axes.
        depth : float
            Length of y axis in nm.
        u : float, optional
            Position on the spline at which local Cartesian coordinates will be built.
        scale: nm, default 1.0
            Scale of coordinates, i.e. spacing of the grid.

        Returns
        -------
        np.ndarray
            (D, V, S, H) shape. Each cooresponds to dimensional vertical, longitudinal
            and horizontal axis, which is ready to be used in `ndi.map_coordinates`.
        """

        mole = self.anchors_to_molecules(u)
        nz = roundint(shape[0] / scale)
        ny = roundint(depth / scale)
        nx = roundint(shape[1] / scale)
        return mole.local_coordinates(shape=(nz, ny, nx), scale=scale)

    def local_cylindrical(
        self,
        r_range: tuple[nm, nm],
        depth: nm,
        u: float | None = None,
        scale: nm = 1.0,
    ) -> NDArray[np.float32]:
        """
        Generate local cylindrical coordinate systems.

        The generated array can be used for `ndi.map_coordinates`. The result coordinate
        systems are flat, i.e., not distorted by the curvature of spline.

        Parameters
        ----------
        r_range : (float, float)
            Lower and upper bound of radius in nm.
        depth : nm
            Length of y axis in nm.
        u : float
            Position on the spline at which local cylindrical coordinates will be built.
        scale: nm, default 1.0
            Scale of coordinates, i.e. spacing of the grid.

        Returns
        -------
        np.ndarray
            (D, V, S, H) shape. Each cooresponds to dimensional, radius, longitudinal
            and angle axis, which is ready to be used in `ndi.map_coordinates`.
        """
        if u is None:
            u = self.anchors
        rmin, rmax = r_range
        ds = self.map(u, der=1)
        ds_norm: NDArray[np.float32] = ds.reshape(-1, 1) / np.sqrt(sum(ds**2))
        depth_px = roundint(depth / scale)
        grid = np.linspace(-depth_px / 2 + 0.5, depth_px / 2 - 0.5, depth_px)
        dy = ds_norm * grid
        y_ax_coords = (self.map(u) / scale).reshape(1, -1) + dy.T
        dslist = np.stack([ds] * depth_px, axis=0)
        map_ = polar_coords_2d(rmin / scale, rmax / scale)
        map_slice = _stack_coords(map_)
        out = _rot_with_vector(map_slice, y_ax_coords, dslist)
        return np.moveaxis(out, -1, 0)

    def cartesian(
        self,
        shape: tuple[nm, nm],
        s_range: tuple[float, float] = (0, 1),
        scale: nm = 1.0,
    ) -> NDArray[np.float32]:
        """
        Generate a Cartesian coordinate system along spline that can be used for
        `ndi.map_coordinate`. Note that this coordinate system is distorted, thus
        does not reflect real geometry (such as distance and derivatives).

        Parameters
        ----------
        shape : (float, float)
            The ZX-shape of output coordinate system. Center of the array will be
            spline curve itself after coodinate transformation.
        s_range : tuple[float, float], default (0, 1)
            Range of spline. Spline coordinate system will be built between
            `spl[s_range[0]]` and `spl[s_range[1]]`.
        scale: nm, default 1.0
            Scale of coordinates, i.e. spacing of the grid.

        Returns
        -------
        np.ndarray
            (V, S, H, D) shape. Each cooresponds to vertical, longitudinal, horizontal
            and dimensional axis.
        """
        dz = roundint(shape[0] / scale)
        dx = roundint(shape[1] / scale)
        return self._get_coords(_cartesian_coords_2d, (dz, dx), s_range, scale)

    def cylindrical(
        self,
        r_range: tuple[nm, nm],
        s_range: tuple[float, float] = (0, 1),
        scale: nm = 1.0,
    ) -> NDArray[np.float32]:
        """
        Generate a cylindrical coordinate system along spline that can be used for
        `ndi.map_coordinate`. Note that this coordinate system is distorted, thus
        does not reflect real geometry (such as distance and derivatives).

        Parameters
        ----------
        r_range : (nm, nm)
            Range of radius in nm. r=0 will be spline curve itself after coodinate
            transformation.
        s_range : tuple[float, float], default (0, 1)
            Range of spline. Spline coordinate system will be built between
            `spl[s_range[0]]` and `spl[s_range[1]]`.
        scale: nm, default 1.0
            Scale of coordinates, i.e. spacing of the grid.

        Returns
        -------
        np.ndarray
            (V, S, H, D) shape. Each cooresponds to radius, longitudinal, angle and
            dimensional axis.
        """
        rmin = r_range[0] / scale
        rmax = r_range[1] / scale
        return self._get_coords(polar_coords_2d, (rmin, rmax), s_range, scale)

    def y_to_position(
        self, y: NDArray[np.float32], nknots: int = 512
    ) -> NDArray[np.float32]:
        """
        Convert y-coordinate to spline position parameter.

        Parameters
        ----------
        y : array-like
            Y coordinates.
        nknots : int, optional
            Number of knots. Increasing the number of knots will increase the accuracy.
        """
        # almost equal to y / self.length()
        u = np.linspace(0, 1, nknots)
        u_tr = _linear_conversion(u, *self._lims)
        dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
        dist = np.concatenate([[0], np.sqrt(dx**2 + dy**2 + dz**2)]).cumsum()
        tck = splrep(dist, u, k=1)
        out = splev(y, tck)
        return out

    def cartesian_to_world(
        self,
        coords: NDArray[np.float32],
        nknots: int = 512,
    ) -> NDArray[np.float32]:
        """
        Inverse Cartesian coordinate mapping, (z', y', x') to world coordinate.

        Parameters
        ----------
        coords : np.ndarray
            Spline Cartesian coordinates. All the coordinates must be in nm unit.

        Returns
        -------
        np.ndarray
            World coordinates.
        """
        ncoords = coords.shape[0]
        _zs = coords[:, 0]
        _us = self.y_to_position(coords[:, 1], nknots=nknots)
        _xs = coords[:, 2]
        _zeros = np.zeros(ncoords, dtype=np.float32)
        coords_ext = np.stack([_zs, _zeros, _xs], axis=1)
        rot = self.get_rotator(_us)
        out = rot.apply(coords_ext) + self.map(_us)

        return out

    def cylindrical_to_world(self, coords: NDArray[np.float32]) -> NDArray[np.float32]:
        """
        Inverse cylindrical coordinate mapping, (r, y, angle) to world coordinate.

        Parameters
        ----------
        coords : np.ndarray
            Cylindrical coordinates. "r" and "y" must be in scale of "nm", while angle
            must be in radian.

        Returns
        -------
        np.ndarray
            World coordinates.
        """
        radius = coords[:, 0]
        y = coords[:, 1]
        theta = coords[:, 2]
        cart_coords = np.stack(
            [radius * np.sin(theta), y, radius * np.cos(theta)], axis=1
        )

        return self.cartesian_to_world(cart_coords)

    def anchors_to_molecules(
        self,
        positions: float | Sequence[float] | None = None,
        rotation: Sequence[float] | None = None,
    ) -> Molecules:
        """
        Convert coordinates of anchors to `Molecules` instance.

        Coordinates of anchors must be in range from 0 to 1. The y-direction of
        `Molecules` always points at the direction of spline and the z- direction always
        in the plane orthogonal to YX-plane.

        Parameters
        ----------
        positions : iterable of float, optional
            Positions. Between 0 and 1. If not given, anchors are used instead.

        Returns
        -------
        Molecules
            Molecules object of points.
        """
        if positions is None:
            positions = self.anchors
        elif np.isscalar(positions):
            positions = [positions]
        pos = self.map(positions)
        yvec = self.map(positions, der=1)
        rot = axes_to_rotator(None, yvec)
        if rotation is not None:
            rotvec = np.zeros((len(rot), 3), dtype=np.float32)
            rotvec[:, 1] = rotation
            rot = rot * Rotation.from_rotvec(rotvec)
        return Molecules(pos=pos, rot=rot, features={Mole.nth: np.arange(len(pos))})

    def cylindrical_to_molecules(
        self,
        coords: NDArray[np.float32],
    ) -> Molecules:
        """
        Convert coordinates of points near the spline to `Molecules` instance.

        Coordinates of points must be those in spline cylindrical coordinate system.

        Parameters
        ----------
        coords : (N, 3) array
            Spline cylindrical coordinates of points.

        Returns
        -------
        Molecules
            Molecules object of points.
        """
        world_coords = self.cylindrical_to_world(coords)

        # world coordinates of the projection point of coords onto the spline
        u = self.y_to_position(coords[:, 1])
        ycoords = self.map(u)
        zvec = world_coords - ycoords
        yvec = self.map(u, der=1)
        return Molecules.from_axes(pos=world_coords, z=zvec, y=yvec)

    def _get_coords(
        self,
        map_func: Callable[[tuple], NDArray[np.float32]],
        map_params: tuple,
        s_range: tuple[float, float],
        scale: nm,
    ):
        """
        Make coordinate system using function `map_func` and stack the same point cloud
        in the direction of the spline, in the range of `s_range`.
        """
        u, y_ax_coords = self._get_y_ax_coords(s_range, scale)
        dslist = self.map(u, 1).astype(np.float32)
        map_ = map_func(*map_params)
        map_slice = _stack_coords(map_)
        out = _rot_with_vector(map_slice, y_ax_coords, dslist)
        return np.moveaxis(out, -1, 0)

    def _get_y_ax_coords(self, s_range: tuple[float, float], scale: nm):
        s0, s1 = s_range
        length = self.length(start=s0, stop=s1)
        stop_length, n_segs = interval_divmod(length, scale)
        n_pixels = n_segs + 1
        s2 = (s1 - s0) * stop_length / length + s0
        if n_pixels < 2:
            raise ValueError("Too short. Change 's_range'.")
        d0, d2 = self.distances([s0, s2])
        u = self.y_to_position(np.linspace(d0, d2, n_pixels))
        y = self.map(u) / scale  # world coordinates of y-axis in spline coords system
        return u, y
anchors deletable property writable

Local anchors along spline.

coeff property

Spline coefficient.

config property

Return the spline configuration

extrapolate property

Extrapolation mode of the spline.

has_anchors property

True if there are any anchors.

knots property

Spline knots.

lims property

Return spline limit positions.

order property

Spline order.

params property

Spline parameters.

props property

Return the spline properties

__repr__()

Use start/end points to describe a spline.

Source code in cylindra/components/spline/_spline_base.py 
282
283
284
285
286
287
def __repr__(self) -> str:
    """Use start/end points to describe a spline."""
    start, end = self.map(self._lims)
    start = "({:.1f}, {:.1f}, {:.1f})".format(*start)
    end = "({:.1f}, {:.1f}, {:.1f})".format(*end)
    return f"Spline[{start}:{end}]"
anchors_to_molecules(positions=None, rotation=None)

Convert coordinates of anchors to Molecules instance.

Coordinates of anchors must be in range from 0 to 1. The y-direction of Molecules always points at the direction of spline and the z- direction always in the plane orthogonal to YX-plane.

Parameters:

Name Type Description Default
positions iterable of float

Positions. Between 0 and 1. If not given, anchors are used instead.

 None

Returns:

Type Description
Molecules

Molecules object of points.
Source code in cylindra/components/spline/_spline_base.py 
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
def anchors_to_molecules(
    self,
    positions: float | Sequence[float] | None = None,
    rotation: Sequence[float] | None = None,
) -> Molecules:
    """
    Convert coordinates of anchors to `Molecules` instance.

    Coordinates of anchors must be in range from 0 to 1. The y-direction of
    `Molecules` always points at the direction of spline and the z- direction always
    in the plane orthogonal to YX-plane.

    Parameters
    ----------
    positions : iterable of float, optional
        Positions. Between 0 and 1. If not given, anchors are used instead.

    Returns
    -------
    Molecules
        Molecules object of points.
    """
    if positions is None:
        positions = self.anchors
    elif np.isscalar(positions):
        positions = [positions]
    pos = self.map(positions)
    yvec = self.map(positions, der=1)
    rot = axes_to_rotator(None, yvec)
    if rotation is not None:
        rotvec = np.zeros((len(rot), 3), dtype=np.float32)
        rotvec[:, 1] = rotation
        rot = rot * Rotation.from_rotvec(rotvec)
    return Molecules(pos=pos, rot=rot, features={Mole.nth: np.arange(len(pos))})
cartesian(shape, s_range=(0, 1), scale=1.0)

Generate a Cartesian coordinate system along spline that can be used for ndi.map_coordinate. Note that this coordinate system is distorted, thus does not reflect real geometry (such as distance and derivatives).

Parameters:

Name Type Description Default
shape (float, float)

The ZX-shape of output coordinate system. Center of the array will be spline curve itself after coodinate transformation.

 required
s_range tuple[float, float]

Range of spline. Spline coordinate system will be built between spl[s_range[0]] and spl[s_range[1]].

 (0, 1)
scale nm

Scale of coordinates, i.e. spacing of the grid.

 1.0

Returns:

Type Description
ndarray

(V, S, H, D) shape. Each cooresponds to vertical, longitudinal, horizontal and dimensional axis.
Source code in cylindra/components/spline/_spline_base.py 
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
def cartesian(
    self,
    shape: tuple[nm, nm],
    s_range: tuple[float, float] = (0, 1),
    scale: nm = 1.0,
) -> NDArray[np.float32]:
    """
    Generate a Cartesian coordinate system along spline that can be used for
    `ndi.map_coordinate`. Note that this coordinate system is distorted, thus
    does not reflect real geometry (such as distance and derivatives).

    Parameters
    ----------
    shape : (float, float)
        The ZX-shape of output coordinate system. Center of the array will be
        spline curve itself after coodinate transformation.
    s_range : tuple[float, float], default (0, 1)
        Range of spline. Spline coordinate system will be built between
        `spl[s_range[0]]` and `spl[s_range[1]]`.
    scale: nm, default 1.0
        Scale of coordinates, i.e. spacing of the grid.

    Returns
    -------
    np.ndarray
        (V, S, H, D) shape. Each cooresponds to vertical, longitudinal, horizontal
        and dimensional axis.
    """
    dz = roundint(shape[0] / scale)
    dx = roundint(shape[1] / scale)
    return self._get_coords(_cartesian_coords_2d, (dz, dx), s_range, scale)
cartesian_to_world(coords, nknots=512)

Inverse Cartesian coordinate mapping, (z', y', x') to world coordinate.

Parameters:

Name Type Description Default
coords ndarray

Spline Cartesian coordinates. All the coordinates must be in nm unit.

 required

Returns:

Type Description
ndarray

World coordinates.
Source code in cylindra/components/spline/_spline_base.py 
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
def cartesian_to_world(
    self,
    coords: NDArray[np.float32],
    nknots: int = 512,
) -> NDArray[np.float32]:
    """
    Inverse Cartesian coordinate mapping, (z', y', x') to world coordinate.

    Parameters
    ----------
    coords : np.ndarray
        Spline Cartesian coordinates. All the coordinates must be in nm unit.

    Returns
    -------
    np.ndarray
        World coordinates.
    """
    ncoords = coords.shape[0]
    _zs = coords[:, 0]
    _us = self.y_to_position(coords[:, 1], nknots=nknots)
    _xs = coords[:, 2]
    _zeros = np.zeros(ncoords, dtype=np.float32)
    coords_ext = np.stack([_zs, _zeros, _xs], axis=1)
    rot = self.get_rotator(_us)
    out = rot.apply(coords_ext) + self.map(_us)

    return out
clip(start, stop)

Clip spline and generate a new one.

This method does not convert spline bases. _lims is updated instead.

Parameters:

Name Type Description Default
start float

New starting position.

 required
stop float

New stopping position.

 required

Returns:

Type Description
Spline

Clipped spline.
Source code in cylindra/components/spline/_spline_base.py 
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
def clip(self, start: float, stop: float) -> Self:
    """
    Clip spline and generate a new one.

    This method does not convert spline bases. `_lims` is updated instead.

    Parameters
    ----------
    start : float
        New starting position.
    stop : float
        New stopping position.

    Returns
    -------
    Spline
        Clipped spline.
    """
    u0 = _linear_conversion(start, *self._lims)
    u1 = _linear_conversion(stop, *self._lims)
    return self.__class__(
        order=self.order,
        lims=(u0, u1),
        extrapolate=self.extrapolate,
        config=self.config,
    )._set_params(self._tck, self._u)
close_to(other)

True if two objects draws the same curve.

Source code in cylindra/components/spline/_spline_base.py 
289
290
291
292
293
294
295
296
297
298
299
300
301
def close_to(self: Self, other: Self) -> bool:
    """True if two objects draws the same curve."""
    if not isinstance(other, self.__class__):
        return False
    t0, c0, k0 = self._tck
    t1, c1, k1 = other._tck
    return (
        np.allclose(t0, t1)
        and all(np.allclose(x, y) for x, y in zip(c0, c1, strict=True))
        and k0 == k1
        and np.allclose(self._u, other._u)
        and np.allclose(self._lims, other._lims)
    )
copy(copy_props=True, copy_config=True)

Copy Spline object.

Parameters:

Name Type Description Default
copy_props bool

Also copy local/global properties if true.

 True

Returns:

Type Description
Spline

Copied object.
Source code in cylindra/components/spline/_spline_base.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def copy(self, copy_props: bool = True, copy_config: bool = True) -> Self:
    """
    Copy Spline object.

    Parameters
    ----------
    copy_props : bool, default True
        Also copy local/global properties if true.

    Returns
    -------
    Spline
        Copied object.
    """
    new = self.__class__(
        order=self.order, lims=self._lims, extrapolate=self.extrapolate
    )
    new._tck = self._tck
    new._u = self._u
    new._anchors = self._anchors

    if copy_props:
        new._props = self.props.copy()
    if copy_config:
        new._config = self.config.copy()
    return new
curvature(positions=None)

Calculate curvature of spline curve.

Parameters:

Name Type Description Default
positions sequence of float

Positions. Between 0 and 1. If not given, anchors are used instead.

 None

Returns:

Type Description
ndarray

Array of curvature.
References
  • https://en.wikipedia.org/wiki/Curvature#Space_curves
Source code in cylindra/components/spline/_spline_base.py 
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
def curvature(
    self,
    positions: Sequence[float] | None = None,
) -> NDArray[np.float32]:
    """
    Calculate curvature of spline curve.

    Parameters
    ----------
    positions : sequence of float, optional
        Positions. Between 0 and 1. If not given, anchors are used instead.

    Returns
    -------
    np.ndarray
        Array of curvature.

    References
    ----------
    - https://en.wikipedia.org/wiki/Curvature#Space_curves
    """

    if positions is None:
        positions = self.anchors

    dz, dy, dx = self.map(positions, der=1).T
    ddz, ddy, ddx = self.map(positions, der=2).T
    a = (
        (ddz * dy - ddy * dz) ** 2
        + (ddx * dz - ddz * dx) ** 2
        + (ddy * dx - ddx * dy) ** 2
    )
    return np.sqrt(a) / (dx**2 + dy**2 + dz**2) ** 1.5
curvature_radii(positions=None)

Inverse of curvature.

Source code in cylindra/components/spline/_spline_base.py 
705
706
707
def curvature_radii(self, positions: Sequence[float] = None) -> NDArray[np.float32]:
    """Inverse of curvature."""
    return 1.0 / self.curvature(positions)
cylindrical(r_range, s_range=(0, 1), scale=1.0)

Generate a cylindrical coordinate system along spline that can be used for ndi.map_coordinate. Note that this coordinate system is distorted, thus does not reflect real geometry (such as distance and derivatives).

Parameters:

Name Type Description Default
r_range (nm, nm)

Range of radius in nm. r=0 will be spline curve itself after coodinate transformation.

 required
s_range tuple[float, float]

Range of spline. Spline coordinate system will be built between spl[s_range[0]] and spl[s_range[1]].

 (0, 1)
scale nm

Scale of coordinates, i.e. spacing of the grid.

 1.0

Returns:

Type Description
ndarray

(V, S, H, D) shape. Each cooresponds to radius, longitudinal, angle and dimensional axis.
Source code in cylindra/components/spline/_spline_base.py 
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
def cylindrical(
    self,
    r_range: tuple[nm, nm],
    s_range: tuple[float, float] = (0, 1),
    scale: nm = 1.0,
) -> NDArray[np.float32]:
    """
    Generate a cylindrical coordinate system along spline that can be used for
    `ndi.map_coordinate`. Note that this coordinate system is distorted, thus
    does not reflect real geometry (such as distance and derivatives).

    Parameters
    ----------
    r_range : (nm, nm)
        Range of radius in nm. r=0 will be spline curve itself after coodinate
        transformation.
    s_range : tuple[float, float], default (0, 1)
        Range of spline. Spline coordinate system will be built between
        `spl[s_range[0]]` and `spl[s_range[1]]`.
    scale: nm, default 1.0
        Scale of coordinates, i.e. spacing of the grid.

    Returns
    -------
    np.ndarray
        (V, S, H, D) shape. Each cooresponds to radius, longitudinal, angle and
        dimensional axis.
    """
    rmin = r_range[0] / scale
    rmax = r_range[1] / scale
    return self._get_coords(polar_coords_2d, (rmin, rmax), s_range, scale)
cylindrical_to_molecules(coords)

Convert coordinates of points near the spline to Molecules instance.

Coordinates of points must be those in spline cylindrical coordinate system.

Parameters:

Name Type Description Default
coords (N, 3) array

Spline cylindrical coordinates of points.

 required

Returns:

Type Description
Molecules

Molecules object of points.
Source code in cylindra/components/spline/_spline_base.py 
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
def cylindrical_to_molecules(
    self,
    coords: NDArray[np.float32],
) -> Molecules:
    """
    Convert coordinates of points near the spline to `Molecules` instance.

    Coordinates of points must be those in spline cylindrical coordinate system.

    Parameters
    ----------
    coords : (N, 3) array
        Spline cylindrical coordinates of points.

    Returns
    -------
    Molecules
        Molecules object of points.
    """
    world_coords = self.cylindrical_to_world(coords)

    # world coordinates of the projection point of coords onto the spline
    u = self.y_to_position(coords[:, 1])
    ycoords = self.map(u)
    zvec = world_coords - ycoords
    yvec = self.map(u, der=1)
    return Molecules.from_axes(pos=world_coords, z=zvec, y=yvec)
cylindrical_to_world(coords)

Inverse cylindrical coordinate mapping, (r, y, angle) to world coordinate.

Parameters:

Name Type Description Default
coords ndarray

Cylindrical coordinates. "r" and "y" must be in scale of "nm", while angle must be in radian.

 required

Returns:

Type Description
ndarray

World coordinates.
Source code in cylindra/components/spline/_spline_base.py 
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
def cylindrical_to_world(self, coords: NDArray[np.float32]) -> NDArray[np.float32]:
    """
    Inverse cylindrical coordinate mapping, (r, y, angle) to world coordinate.

    Parameters
    ----------
    coords : np.ndarray
        Cylindrical coordinates. "r" and "y" must be in scale of "nm", while angle
        must be in radian.

    Returns
    -------
    np.ndarray
        World coordinates.
    """
    radius = coords[:, 0]
    y = coords[:, 1]
    theta = coords[:, 2]
    cart_coords = np.stack(
        [radius * np.sin(theta), y, radius * np.cos(theta)], axis=1
    )

    return self.cartesian_to_world(cart_coords)
distances(positions=None, nknots=512)

Get the distances from u=0.

Parameters:

Name Type Description Default
positions sequence of float

Positions. Between 0 and 1. If not given, anchors are used instead.

 None

Returns:

Type Description
ndarray

Distances for each u.
Source code in cylindra/components/spline/_spline_base.py 
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
def distances(
    self,
    positions: Sequence[float] | None = None,
    nknots: int = 512,
) -> NDArray[np.float32]:
    """
    Get the distances from u=0.

    Parameters
    ----------
    positions : sequence of float, optional
        Positions. Between 0 and 1. If not given, anchors are used instead.

    Returns
    -------
    np.ndarray
        Distances for each `u`.
    """
    if positions is None:
        _u = self.anchors
    else:
        _u = np.asarray(positions, dtype=np.float32)
        if _u.ndim != 1:
            raise ValueError("Positions must be 1D array.")
    u = np.linspace(0, 1, nknots)
    u_tr = _linear_conversion(u, *self._lims)
    dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
    dist = np.concatenate([[0], np.sqrt(dx**2 + dy**2 + dz**2)]).cumsum()
    tck = splrep(u, dist, k=1)
    out = splev(_u, tck)
    return out
fit(coords, *, err_max=1.0)

Fit spline model to coordinates.

This method uses scipy.interpolate.splprep to fit given coordinates to a spline.

Parameters:

Name Type Description Default
coords ndarray

Coordinates. Must be (N, 3).

 required
err_max float

Error allowed for fitting. Several upper limit of residual values will be used and the fit that results in error lower than this value and minimize the maximum curvature will be chosen.

 1.0

Returns:

Type Description
Spline

New spline fit to given coordinates.
Source code in cylindra/components/spline/_spline_base.py 
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
def fit(
    self,
    coords: ArrayLike,
    *,
    err_max: nm = 1.0,
) -> Self:
    """
    Fit spline model to coordinates.

    This method uses `scipy.interpolate.splprep` to fit given coordinates to a
    spline.

    Parameters
    ----------
    coords : np.ndarray
        Coordinates. Must be (N, 3).
    err_max : float, default 1.0
        Error allowed for fitting. Several upper limit of residual values will be
        used and the fit that results in error lower than this value and minimize
        the maximum curvature will be chosen.

    Returns
    -------
    Spline
        New spline fit to given coordinates.
    """
    crds = np.asarray(coords)
    npoints = crds.shape[0]
    if npoints < 2:
        raise ValueError("Number of input coordinates must be > 1.")
    if npoints <= self.order:
        k = npoints - 1
    else:
        k = self.order

    if self.is_inverted():
        crds = crds[::-1]

    if err_max > 4.0:
        raise ValueError("std_max must be smaller than 4.0.")
    if err_max < 1e-3:
        std_list = [err_max]
    else:
        ntrial = max(int(err_max / 0.02), 2)
        std_list = np.linspace(0, err_max, ntrial)[1:]

    fit_results = list[SplineFitResult]()
    new = self.__class__(order=k, extrapolate=self.extrapolate, config=self.config)
    with warnings.catch_warnings():
        # fitting may fail for some std
        warnings.simplefilter("ignore", RuntimeWarning)
        for std in std_list:
            _tck, _u = splprep(crds.T, k=k, s=std**2 * npoints)
            new._set_params(_tck, _u)
            _crds_at_u = new.map(_u)
            res: NDArray[np.float32] = np.sqrt(
                np.sum((_crds_at_u - crds) ** 2, axis=1)
            )
            _knots = _tck[0][new.order : -new.order]
            nedge = _knots.size - 1
            assert nedge > 0
            nanc = nedge * 20 + 1
            anc = np.interp(
                np.linspace(0, 1, nanc), np.linspace(0, 1, nedge + 1), _knots
            )
            max_curvature = new.curvature(anc).max()
            success = res.max() <= err_max
            fit_results.append(
                SplineFitResult((_tck, _u), max_curvature, res, success)
            )

    fit_results_filt = list(filter(lambda x: x.success, fit_results))
    if len(fit_results_filt) == 0:
        fit_results_filt = fit_results

    reult_opt = min(fit_results_filt, key=lambda x: x.curvature)
    return new._set_params(*reult_opt.params)
from_dict(d) classmethod

Construct a spline model from a dictionary.

Parameters:

Name Type Description Default
d SplineInfo

Dictionary with keys "t", "c", "k", "u" and "lims".

 required

Returns:

Type Description
Spline

Spline object constructed from the dictionary.
Source code in cylindra/components/spline/_spline_base.py 
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
@classmethod
def from_dict(cls: type[Self], d: SplineInfo) -> Self:
    """
    Construct a spline model from a dictionary.

    Parameters
    ----------
    d: dict
        Dictionary with keys "t", "c", "k", "u" and "lims".

    Returns
    -------
    Spline
        Spline object constructed from the dictionary.
    """
    self = cls(
        order=d.get("k", 3),
        lims=d.get("lims", (0, 1)),
        extrapolate=d.get("extrapolate", "linear"),
    )
    t = np.asarray(d["t"])
    c = [np.asarray(d["c"][k]) for k in "zyx"]
    k = roundint(d["k"])
    self._tck = (t, c, k)
    self._u = np.asarray(d["u"])
    self.props._window_size = d.get("localprops_window_size", {})
    self.props._binsize_loc = d.get("binsize_loc", {})
    self.props._binsize_glob = d.get("binsize_glob", {})
    if cfg := d.get("config", None):
        self._config = SplineConfig.from_dict(cfg)
    return self
get_rotator(positions=None, inverse=False)

Calculate list of Affine transformation matrix along spline, which correspond to the orientation of spline curve.

Parameters:

Name Type Description Default
positions (array - like, (N,))

Positions. Between 0 and 1.

 None
inverse bool

If True, rotation matrix will be inversed.

 False

Returns:

Type Description
Rotation

Rotation object at each anchor.
Source code in cylindra/components/spline/_spline_base.py 
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
def get_rotator(
    self,
    positions: Sequence[float] | None = None,
    inverse: bool = False,
) -> Rotation:
    """
    Calculate list of Affine transformation matrix along spline, which correspond to
    the orientation of spline curve.

    Parameters
    ----------
    positions : array-like, (N,)
        Positions. Between 0 and 1.
    inverse : bool, default False
        If True, rotation matrix will be inversed.

    Returns
    -------
    Rotation
        Rotation object at each anchor.
    """
    if positions is None:
        positions = self.anchors
    ds = self.map(positions, der=1)
    out = axes_to_rotator(None, -ds)

    if inverse:
        out = out.inv()

    return out
has_props()

True if there are any properties.

Source code in cylindra/components/spline/_spline_base.py 
81
82
83
def has_props(self) -> bool:
    """True if there are any properties."""
    return len(self.props.loc) > 0 or len(self.props.glob) > 0
invert()

Invert the direction of spline.

Returns:

Type Description
Spline

Inverted object
Source code in cylindra/components/spline/_spline_base.py 
627
628
629
630
631
632
633
634
635
636
637
638
639
640
def invert(self) -> Self:
    """
    Invert the direction of spline.

    Returns
    -------
    Spline
        Inverted object
    """
    anchors = self._anchors
    inverted = self.clip(1.0, 0.0)
    if anchors is not None:
        inverted.anchors = 1 - anchors[::-1]
    return inverted
is_inverted()

Return true if spline is inverted.

Source code in cylindra/components/spline/_spline_base.py 
224
225
226
def is_inverted(self) -> bool:
    """Return true if spline is inverted."""
    return self._lims[0] > self._lims[1]
length(start=0, stop=1, nknots=512)

Approximate the length of B-spline between [start, stop] by partitioning the spline with 'nknots' knots. nknots=256 is large enough for most cases.

Source code in cylindra/components/spline/_spline_base.py 
616
617
618
619
620
621
622
623
624
625
def length(self, start: float = 0, stop: float = 1, nknots: int = 512) -> nm:
    """
    Approximate the length of B-spline between [start, stop] by partitioning
    the spline with 'nknots' knots. nknots=256 is large enough for most cases.
    """
    _assert_fitted(self)
    u = np.linspace(start, stop, nknots)
    u_tr = _linear_conversion(u, *self._lims)
    dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
    return np.sum(np.sqrt(dx**2 + dy**2 + dz**2))
line(start, end) classmethod

Create a line spline.

Parameters:

Name Type Description Default
start array_like

Start point of the line.

 required
end array_like

End point of the line.

 required

Returns:

Type Description
Spline

Line spline.
Source code in cylindra/components/spline/_spline_base.py 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
@classmethod
def line(cls, start: ArrayLike, end: ArrayLike) -> Self:
    """
    Create a line spline.

    Parameters
    ----------
    start : array_like
        Start point of the line.
    end : array_like
        End point of the line.

    Returns
    -------
    Spline
        Line spline.
    """
    spl = cls()
    coords = np.stack([start, end], axis=0)
    return spl.fit(coords, err_max=0.0)
local_cartesian(shape, depth, u=None, scale=1.0)

Generate local Cartesian coordinate systems.

The generated array can be used for ndi.map_coordinates. The result coordinate systems are flat, i.e., not distorted by the curvature of spline.

Parameters:

Name Type Description Default
shape (float, float)

Vertical and horizontal length of Cartesian coordinates. Corresponds to zx axes.

 required
depth float

Length of y axis in nm.

 required
u float

Position on the spline at which local Cartesian coordinates will be built.

 None
scale nm

Scale of coordinates, i.e. spacing of the grid.

 1.0

Returns:

Type Description
ndarray

(D, V, S, H) shape. Each cooresponds to dimensional vertical, longitudinal and horizontal axis, which is ready to be used in ndi.map_coordinates.
Source code in cylindra/components/spline/_spline_base.py 
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
def local_cartesian(
    self,
    shape: tuple[nm, nm],
    depth: nm,
    u: float | Sequence[float] = None,
    scale: nm = 1.0,
) -> NDArray[np.float32]:
    """
    Generate local Cartesian coordinate systems.

    The generated array can be used for `ndi.map_coordinates`. The result coordinate
    systems are flat, i.e., not distorted by the curvature of spline.

    Parameters
    ----------
    shape : (float, float)
        Vertical and horizontal length of Cartesian coordinates. Corresponds to zx
        axes.
    depth : float
        Length of y axis in nm.
    u : float, optional
        Position on the spline at which local Cartesian coordinates will be built.
    scale: nm, default 1.0
        Scale of coordinates, i.e. spacing of the grid.

    Returns
    -------
    np.ndarray
        (D, V, S, H) shape. Each cooresponds to dimensional vertical, longitudinal
        and horizontal axis, which is ready to be used in `ndi.map_coordinates`.
    """

    mole = self.anchors_to_molecules(u)
    nz = roundint(shape[0] / scale)
    ny = roundint(depth / scale)
    nx = roundint(shape[1] / scale)
    return mole.local_coordinates(shape=(nz, ny, nx), scale=scale)
local_cylindrical(r_range, depth, u=None, scale=1.0)

Generate local cylindrical coordinate systems.

The generated array can be used for ndi.map_coordinates. The result coordinate systems are flat, i.e., not distorted by the curvature of spline.

Parameters:

Name Type Description Default
r_range (float, float)

Lower and upper bound of radius in nm.

 required
depth nm

Length of y axis in nm.

 required
u float

Position on the spline at which local cylindrical coordinates will be built.

 None
scale nm

Scale of coordinates, i.e. spacing of the grid.

 1.0

Returns:

Type Description
ndarray

(D, V, S, H) shape. Each cooresponds to dimensional, radius, longitudinal and angle axis, which is ready to be used in ndi.map_coordinates.
Source code in cylindra/components/spline/_spline_base.py 
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
def local_cylindrical(
    self,
    r_range: tuple[nm, nm],
    depth: nm,
    u: float | None = None,
    scale: nm = 1.0,
) -> NDArray[np.float32]:
    """
    Generate local cylindrical coordinate systems.

    The generated array can be used for `ndi.map_coordinates`. The result coordinate
    systems are flat, i.e., not distorted by the curvature of spline.

    Parameters
    ----------
    r_range : (float, float)
        Lower and upper bound of radius in nm.
    depth : nm
        Length of y axis in nm.
    u : float
        Position on the spline at which local cylindrical coordinates will be built.
    scale: nm, default 1.0
        Scale of coordinates, i.e. spacing of the grid.

    Returns
    -------
    np.ndarray
        (D, V, S, H) shape. Each cooresponds to dimensional, radius, longitudinal
        and angle axis, which is ready to be used in `ndi.map_coordinates`.
    """
    if u is None:
        u = self.anchors
    rmin, rmax = r_range
    ds = self.map(u, der=1)
    ds_norm: NDArray[np.float32] = ds.reshape(-1, 1) / np.sqrt(sum(ds**2))
    depth_px = roundint(depth / scale)
    grid = np.linspace(-depth_px / 2 + 0.5, depth_px / 2 - 0.5, depth_px)
    dy = ds_norm * grid
    y_ax_coords = (self.map(u) / scale).reshape(1, -1) + dy.T
    dslist = np.stack([ds] * depth_px, axis=0)
    map_ = polar_coords_2d(rmin / scale, rmax / scale)
    map_slice = _stack_coords(map_)
    out = _rot_with_vector(map_slice, y_ax_coords, dslist)
    return np.moveaxis(out, -1, 0)
make_anchors(interval=None, n=None, max_interval=None)

Make anchor points at constant intervals. Either interval, number of anchor or the maximum interval between anchors can be specified.

Parameters:

Name Type Description Default
interval nm

Interval between anchor points.

 None
n int

Number of anchor points, including both ends.

 None
max_interval nm | None

Spline will be split by as little anchors as possible but interval between anchors will not be larger than this. The number of anchors are also guaranteed to be larger than spline order.

 None
Source code in cylindra/components/spline/_spline_base.py 
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def make_anchors(
    self,
    interval: nm | None = None,
    n: int | None = None,
    max_interval: nm | None = None,
) -> Self:
    """
    Make anchor points at constant intervals. Either interval, number of anchor or
    the maximum interval between anchors can be specified.

    Parameters
    ----------
    interval : nm, optional
        Interval between anchor points.
    n : int, optional
        Number of anchor points, including both ends.
    max_interval: nm, optional
        Spline will be split by as little anchors as possible but interval between
        anchors will not be larger than this. The number of anchors are also
        guaranteed to be larger than spline order.
    """
    self.anchors = self.prep_anchor_positions(interval, n, max_interval)
    return self
map(positions=None, der=0)

Calculate coordinates (or n-th derivative) at points on the spline.

Parameters:

Name Type Description Default
positions ndarray or float

Positions. Between 0 and 1. If not given, anchors are used instead.

 None
der int

der-th derivative will be calculated.

 0

Returns:

Type Description
ndarray

Positions or vectors in (3,) or (N, 3) shape.
Source code in cylindra/components/spline/_spline_base.py 
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
def map(
    self,
    positions: float | NDArray[np.number] | None = None,
    der: int = 0,
) -> NDArray[np.float32]:
    """
    Calculate coordinates (or n-th derivative) at points on the spline.

    Parameters
    ----------
    positions : np.ndarray or float, optional
        Positions. Between 0 and 1. If not given, anchors are used instead.
    der : int, default 0
        `der`-th derivative will be calculated.

    Returns
    -------
    np.ndarray
        Positions or vectors in (3,) or (N, 3) shape.
    """
    _assert_fitted(self)
    if positions is None:
        positions = self.anchors
    u0, u1 = self._lims
    if np.isscalar(positions):
        if self.order < der:
            return np.zeros(3, dtype=np.float32)
        u_tr = _linear_conversion(float(positions), u0, u1)
        if 0 <= u_tr <= 1 or self.extrapolate is ExtrapolationMode.default:
            coord = splev([u_tr], self._tck, der=der)
        elif self.extrapolate is ExtrapolationMode.linear:
            if der == 0:
                if u_tr < 0:
                    der0 = splev([0], self._tck, der=0)
                    der1 = splev([0], self._tck, der=1)
                    dr = u_tr
                else:
                    der0 = splev([1], self._tck, der=0)
                    der1 = splev([1], self._tck, der=1)
                    dr = u_tr - 1
                coord = [a0 + a1 * dr for a0, a1 in zip(der0, der1, strict=True)]
            elif der == 1:
                if u_tr < 0:
                    coord = splev([0], self._tck, der=1)
                else:
                    coord = splev([1], self._tck, der=1)
            else:
                coord = [[0], [0], [0]]
        else:
            raise ValueError(f"Invalid extrapolation mode: {self.extrapolate!r}.")
        out = np.concatenate(coord).astype(np.float32)

    else:
        u_tr = _linear_conversion(np.asarray(positions, dtype=np.float32), u0, u1)
        if self.order < der:
            return np.zeros((u_tr.size, 3), dtype=np.float32)
        if self.extrapolate is ExtrapolationMode.default:
            out = np.stack(splev(u_tr, self._tck, der=der), axis=1).astype(
                np.float32
            )
        elif self.extrapolate is ExtrapolationMode.linear:
            sl_small = u_tr < 0
            sl_large = u_tr > 1
            n_small = np.count_nonzero(sl_small)
            n_large = np.count_nonzero(sl_large)
            out = np.stack(splev(u_tr, self._tck, der=der), axis=1).astype(
                np.float32
            )
            if der == 0:
                if n_small > 0:
                    der0 = np.array(splev(0, self._tck, der=0), dtype=np.float32)
                    der1 = np.array(splev(0, self._tck, der=1), dtype=np.float32)
                    dr = u_tr[sl_small]
                    coords_new = der0 + der1 * dr[:, np.newaxis]
                    out[sl_small] = coords_new

                if n_large > 0:
                    der0 = splev(1, self._tck, der=0)
                    der1 = splev(1, self._tck, der=1)
                    dr = u_tr[sl_large] - 1
                    coords_new = der0 + der1 * dr[:, np.newaxis]
                    out[sl_large] = coords_new

            elif der == 1:
                if n_small > 0:
                    out[sl_small] = splev(0, self._tck, der=1)
                if n_large > 0:
                    out[sl_large] = splev(1, self._tck, der=1)
            else:
                if n_small > 0:
                    out[sl_small] = 0
                if n_large > 0:
                    out[sl_large] = 0
        else:
            raise ValueError(f"Invalid extrapolation mode: {self.extrapolate!r}.")

    if u0 > u1 and der % 2 == 1:
        out = -out
    return out
partition(n, der=0)

Return the n-partitioning coordinates of the spline.

Source code in cylindra/components/spline/_spline_base.py 
611
612
613
614
def partition(self, n: int, der: int = 0) -> NDArray[np.float32]:
    """Return the n-partitioning coordinates of the spline."""
    u = np.linspace(0, 1, n)
    return self.map(u, der)
resample(max_interval=1.0, err_max=0.1)

Resample a new spline along the original spline.

Parameters:

Name Type Description Default
max_interval nm

Maximum interval between resampling points.

 1.0
err_max float

Spline fitting maximum error.

 0.1

Returns:

Type Description
Spline

Resampled spline object.
Source code in cylindra/components/spline/_spline_base.py 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def resample(self, max_interval: nm = 1.0, err_max: nm = 0.1) -> Self:
    """
    Resample a new spline along the original spline.

    Parameters
    ----------
    max_interval : nm, default 1.0
        Maximum interval between resampling points.
    err_max : float, default 0.1
        Spline fitting maximum error.

    Returns
    -------
    Spline
        Resampled spline object.
    """
    l = self.length()
    points = self.map(np.linspace(0, 1, ceilint(l / max_interval)))
    return self.fit(points, err_max=err_max)
restore()

Restore the original, not-clipped spline.

Returns:

Type Description
Spline

Copy of the original spline.
Source code in cylindra/components/spline/_spline_base.py 
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
def restore(self) -> Self:
    """
    Restore the original, not-clipped spline.

    Returns
    -------
    Spline
        Copy of the original spline.
    """
    return self.__class__(
        order=self.order,
        lims=(0, 1),
        extrapolate=self.extrapolate,
        config=self.config,
    )._set_params(self._tck, self._u)
shift(positions=None, shifts=None, *, err_max=1.0)

Fit spline model using a list of shifts in XZ-plane.

Parameters:

Name Type Description Default
positions sequence of float

Positions. Between 0 and 1. If not given, anchors are used instead.

 None
shifts ndarray

Shift from center in nm. Must be (N, 2).

 None
err_max float

Error allowed for fitting. See Spline.fit.

 1.0

Returns:

Type Description
Spline

Spline shifted by fitting to given coordinates.
Source code in cylindra/components/spline/_spline_base.py 
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
def shift(
    self,
    positions: Sequence[float] | None = None,
    shifts: NDArray[np.floating] | None = None,
    *,
    err_max: nm = 1.0,
) -> Self:
    """
    Fit spline model using a list of shifts in XZ-plane.

    Parameters
    ----------
    positions : sequence of float, optional
        Positions. Between 0 and 1. If not given, anchors are used instead.
    shifts : np.ndarray
        Shift from center in nm. Must be (N, 2).
    err_max : float, default 1.0
        Error allowed for fitting. See `Spline.fit`.

    Returns
    -------
    Spline
        Spline shifted by fitting to given coordinates.
    """
    if shifts is None:
        raise ValueError("Shifts must be given.")
    coords = self.map(positions)
    rot = self.get_rotator(positions)
    # insert 0 in y coordinates.
    shifts = np.stack([shifts[:, 0], np.zeros(len(rot)), shifts[:, 1]], axis=1)
    coords += rot.apply(shifts)
    return self.fit(coords, err_max=err_max)
to_dict()

Convert spline info into a dict.

Source code in cylindra/components/spline/_spline_base.py 
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
def to_dict(self) -> SplineInfo:
    """Convert spline info into a dict."""
    t, c, k = self._tck
    u = self._u
    return {
        "t": t.tolist(),
        "c": {"z": c[0].tolist(), "y": c[1].tolist(), "x": c[2].tolist()},
        "k": k,
        "u": u.tolist(),
        "lims": self._lims,
        "localprops_window_size": dict(self.props.window_size),
        "binsize_loc": dict(self.props.binsize_loc),
        "binsize_glob": dict(self.props.binsize_glob),
        "extrapolate": self._extrapolate.name,
        "config": self.config.asdict(),
    }
translate(shift)

Translate the spline by given shift vectors.

Source code in cylindra/components/spline/_spline_base.py 
178
179
180
181
182
183
def translate(self, shift: tuple[nm, nm, nm]):
    """Translate the spline by given shift vectors."""
    new = self.copy()
    c = [x + s for x, s in zip(self.coeff, shift, strict=True)]
    new._tck = (self.knots, c, self.order)
    return new
with_config(config, copy_props=False)

Return a copy of the spline with a new config.

Source code in cylindra/components/spline/_spline_base.py 
120
121
122
123
124
125
126
127
128
129
130
def with_config(
    self,
    config: dict[str, Any] | SplineConfig,
    copy_props: bool = False,
) -> Self:
    """Return a copy of the spline with a new config."""
    new = self.copy(copy_props=copy_props, copy_config=False)
    if not isinstance(config, SplineConfig):
        config = SplineConfig.construct(**config)
    new._config = config
    return new
with_extrapolation(extrapolate)

Return a copy of the spline with a new extrapolation mode.

Source code in cylindra/components/spline/_spline_base.py 
114
115
116
117
118
def with_extrapolation(self, extrapolate: ExtrapolationMode | str) -> Self:
    """Return a copy of the spline with a new extrapolation mode."""
    new = self.copy()
    new._extrapolate = ExtrapolationMode(extrapolate)
    return new
y_to_position(y, nknots=512)

Convert y-coordinate to spline position parameter.

Parameters:

Name Type Description Default
y array - like

Y coordinates.

 required
nknots int

Number of knots. Increasing the number of knots will increase the accuracy.

 512
Source code in cylindra/components/spline/_spline_base.py 
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
def y_to_position(
    self, y: NDArray[np.float32], nknots: int = 512
) -> NDArray[np.float32]:
    """
    Convert y-coordinate to spline position parameter.

    Parameters
    ----------
    y : array-like
        Y coordinates.
    nknots : int, optional
        Number of knots. Increasing the number of knots will increase the accuracy.
    """
    # almost equal to y / self.length()
    u = np.linspace(0, 1, nknots)
    u_tr = _linear_conversion(u, *self._lims)
    dz, dy, dx = map(np.diff, splev(u_tr, self._tck, der=0))
    dist = np.concatenate([[0], np.sqrt(dx**2 + dy**2 + dz**2)]).cumsum()
    tck = splrep(dist, u, k=1)
    out = splev(y, tck)
    return out

SplineConfig dataclass

Class for spline configuration.

Source code in cylindra/components/spline/_config.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
@dataclass(frozen=True)
class SplineConfig:
    """Class for spline configuration."""

    npf_range: Range[int] = Range(11, 17)
    spacing_range: Range[nm] = Range(3.9, 4.3)
    twist_range: Range[float] = Range(-1.0, 1.0)
    rise_range: Range[float] = Range(0.0, 45.0)
    rise_sign: Literal[-1, 1] = -1
    clockwise: Literal["PlusToMinus", "MinusToPlus"] = "MinusToPlus"
    thickness_inner: nm = 2.0
    thickness_outer: nm = 3.0
    fit_depth: nm = 48.0
    fit_width: nm = 44.0

    def _repr_pretty_(self, p, cycle: bool):
        if cycle:
            p.text(repr(self))
        parts = list[str]()
        for k in self.__dataclass_fields__:
            v = getattr(self, k)
            parts.append(f"{k}={v!r}")
        cont = ",\n\t".join(parts)
        p.text(f"SplineConfig(\n\t{cont}\n)")

    def copy(self) -> SplineConfig:
        return SplineConfig(
            npf_range=self.npf_range.copy(),
            spacing_range=self.spacing_range.copy(),
            twist_range=self.twist_range.copy(),
            rise_range=self.rise_range.copy(),
            rise_sign=self.rise_sign,
            clockwise=self.clockwise,
            thickness_inner=self.thickness_inner,
            thickness_outer=self.thickness_outer,
            fit_depth=self.fit_depth,
            fit_width=self.fit_width,
        )

    def asdict(self) -> dict[str, Any]:
        return {
            "npf_range": self.npf_range.astuple(),
            "spacing_range": self.spacing_range.astuple_rounded(),
            "twist_range": self.twist_range.astuple_rounded(),
            "rise_range": self.rise_range.astuple_rounded(),
            "rise_sign": self.rise_sign,
            "clockwise": self.clockwise,
            "thickness_inner": self.thickness_inner,
            "thickness_outer": self.thickness_outer,
            "fit_depth": self.fit_depth,
            "fit_width": self.fit_width,
        }

    def json_dumps(self) -> str:
        strings = list[str]()
        for k, v in self.asdict().items():
            if isinstance(v, tuple):
                strings.append(f'"{k}": {list(v)!r}'.replace("'", '"'))
            else:
                strings.append(f'"{k}": {v!r}'.replace("'", '"'))
        return "{\n    " + ",\n    ".join(strings) + "\n}"

    @classmethod
    def construct(cls, **kwargs) -> SplineConfig:
        """Construct a SplineConfig with argument check."""
        return SplineConfig().updated(**kwargs)

    @classmethod
    def from_dict(
        cls,
        cfg: dict[str, Any],
        unknown: Literal["warn", "error", "ignore"] = "warn",
    ) -> SplineConfig:
        # for version compatibility
        _undef = {}
        cfg_input = {}
        for k, v in cfg.items():
            if k not in SplineConfig.__annotations__:
                _undef[k] = v
            else:
                cfg_input[k] = v

        if _undef:
            msg = f"Unknown keys, maybe due to version incompatibility: {_undef!r}"
            match unknown:
                case "error":
                    raise ValueError(msg)
                case "warn":
                    warnings.warn(msg, RuntimeWarning, stacklevel=2)
                case "ignore":
                    pass
                case other:  # pragma: no cover
                    raise ValueError(f"Got invalid case {other!r}")
        return SplineConfig().updated(**cfg_input)

    @classmethod
    def from_file(
        cls,
        path: str,
        unknown: Literal["warn", "error", "ignore"] = "warn",
    ) -> SplineConfig:
        with open(path) as f:
            cfg: dict = json.load(f)
        return SplineConfig.from_dict(cfg, unknown=unknown)

    def to_file(self, path: str):
        with open(path, mode="w") as f:
            json.dump(self.asdict(), f)
        return None

    def updated(
        self,
        npf_range: tuple[int, int] | None = None,
        spacing_range: tuple[nm, nm] | None = None,
        twist_range: tuple[float, float] | None = None,
        rise_range: tuple[float, float] | None = None,
        rise_sign: Literal[-1, 1] | None = None,
        clockwise: Literal["PlusToMinus", "MinusToPlus"] | None = None,
        thickness_inner: nm | None = None,
        thickness_outer: nm | None = None,
        fit_depth: nm | None = None,
        fit_width: nm | None = None,
    ) -> SplineConfig:
        kwargs = locals()
        kwargs.pop("self")
        for k, v in kwargs.items():
            if v is None:
                kwargs[k] = getattr(self, k)
        for rng in ["npf_range", "spacing_range", "twist_range", "rise_range"]:
            kwargs[rng] = _norm_range(kwargs[rng])
        if kwargs["rise_sign"] not in [-1, 1]:
            raise ValueError("rise_sign must be -1 or 1")
        if kwargs["clockwise"] not in ["PlusToMinus", "MinusToPlus"]:
            raise ValueError("clockwise must be PlusToMinus or MinusToPlus")
        for n in ["thickness_inner", "thickness_outer", "fit_depth", "fit_width"]:
            if kwargs[n] < 0:
                raise ValueError(f"{n} must be non-negative")
        return SplineConfig(**kwargs)
construct(**kwargs) classmethod

Construct a SplineConfig with argument check.

Source code in cylindra/components/spline/_config.py 
120
121
122
123
@classmethod
def construct(cls, **kwargs) -> SplineConfig:
    """Construct a SplineConfig with argument check."""
    return SplineConfig().updated(**kwargs)

CylTomogram

Bases: Tomogram

Tomogram with cylindrical splines.

Source code in cylindra/components/tomogram/_cyl_tomo.py 
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
class CylTomogram(Tomogram):
    """Tomogram with cylindrical splines."""

    def __init__(self):
        super().__init__()
        self._splines = SplineList()

    @property
    def splines(self) -> SplineList:
        """List of splines."""
        return self._splines

    def add_spline(
        self,
        coords: ArrayLike,
        *,
        order: int = 3,
        err_max: nm = 0.5,
        extrapolate: ExtrapolationMode | str = ExtrapolationMode.linear,
        config: SplineConfig | dict[str, Any] = {},
    ) -> None:
        """
        Add spline path to tomogram.

        Parameters
        ----------
        coords : array-like
            (N, 3) array of coordinates. A spline curve that fit it well is added.
        order : int, optional
            Order of spline curve.
        extrapolate : str, optional
            Extrapolation mode of the spline.
        config : SplineConfig or dict, optional
            Configuration for spline fitting.
        """
        _coords = np.asarray(coords)
        ncoords = _coords.shape[0]
        spl = CylSpline(
            order=order,
            config=config,
            extrapolate=extrapolate,
        ).fit(_coords, err_max=err_max)
        interval: nm = 30.0
        length = spl.length()

        n = int(length / interval) + 1
        fit = spl.map(np.linspace(0, 1, n))
        if ncoords <= spl.order and ncoords < fit.shape[0]:
            return self.add_spline(
                fit,
                order=order,
                err_max=err_max,
                extrapolate=extrapolate,
                config=config,
            )

        self.splines.append(spl)
        return None

    @_misc.batch_process
    def make_anchors(
        self,
        i: int = None,
        *,
        interval: nm | None = None,
        n: int | None = None,
        max_interval: nm | None = None,
    ):
        """
        Make anchors on spline object(s).

        Parameters
        ----------
        interval : nm, optional
            Anchor intervals.
        n : int, optional
            Number of anchors
        max_interval : nm, optional
            Maximum interval between anchors.

        """
        self.splines[i].make_anchors(interval=interval, n=n, max_interval=max_interval)
        return None

    def align_to_polarity(self, orientation: Ori | str = Ori.MinusToPlus) -> Self:
        """
        Align all the splines in the direction parallel to the given polarity.

        Parameters
        ----------
        orientation : Ori or str, default Ori.MinusToPlus
            To which direction splines will be aligned.

        Returns
        -------
        Tomogram object
            Same object with updated splines.
        """
        orientation = Ori(orientation)
        if orientation is Ori.none:
            raise ValueError("Must be PlusToMinus or MinusToPlus.")
        for i, spl in enumerate(self.splines):
            if spl.orientation is Ori.none:
                raise ValueError(f"Spline-{i} has no orientation.")
            if spl.orientation != orientation:
                try:
                    self.splines[i] = spl.invert()
                except Exception as e:
                    raise type(e)(f"Cannot invert spline-{i}: {e}")
        return self

    @_misc.batch_process
    def fit(
        self,
        i: int = None,
        *,
        max_interval: nm = 30.0,
        degree_precision: float = 0.5,
        binsize: int = 1,
        err_max: nm = 1.0,
        edge_sigma: nm = 2.0,
        max_shift: nm = 5.0,
        n_rotations: int = 5,
    ) -> _misc.FitResult:
        """
        Roughly fit splines to cylindrical structures.

        Subtomograms will be sampled at every ``max_interval`` nm. In dense mode,
        Subtomograms will be masked relative to XY-plane, using sigmoid function.
        Sharpness of the sigmoid function is determined by ``dense_mode_sigma``
        (``dense_mode_sigma=0`` corresponds to a step function).

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to fit.
        max_interval : nm, default 30.0
            Maximum interval of sampling points in nm unit.
        degree_precision : float, default 0.5
            Precision of xy-tilt degree in angular correlation.
        binsize : int, default 1
            Multiscale bin size used for fitting.
        edge_sigma : nm, default 2.0
            Sharpness of mask at the edges. If not None, fitting will be executed after regions
            outside the cylinder are masked. Soft mask is important for precision because sharp
            changes in intensity cause strong correlation at the edges.
        max_shift: nm, default 5.0
            Maximum shift from the true center of the cylinder. This parameter is used in phase
            cross correlation.
        n_rotations : int, default 5
            Number of rotations to be tested during finding the cylinder center.

        Returns
        -------
        FitResult
            Result of fitting.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.fit, i={i}")
        spl = self.splines[i]
        anc = spl.prep_anchor_positions(max_interval=max_interval)
        subtomograms, interval, scale = self._prep_fit_spline(spl, anc, binsize)

        with set_gpu():
            subtomograms = _misc.soft_mask_edges(
                subtomograms, spl, anc, scale, edge_sigma
            )
            ds = spl.map(anc, der=1)
            yx_tilt = np.rad2deg(np.arctan2(-ds[:, 2], ds[:, 1]))
            degree_max = 14.0
            nrots = roundint(degree_max / degree_precision) + 1

            # Angular correlation
            out = _misc.dask_angle_corr(subtomograms, yx_tilt, nrots=nrots)
            refined_tilt_deg = np.array(out)
            refined_tilt_rad = np.deg2rad(refined_tilt_deg)

            # If subtomograms are sampled at short intervals, angles should be smoothened to
            # avoid overfitting.
            size = 2 * roundint(48.0 / interval) + 1
            if size > 1:
                # Mirror-mode padding is "a b c d | c b a".
                refined_tilt_rad = _misc.angle_uniform_filter(
                    refined_tilt_rad, size=size, mode=Mode.mirror
                )
                refined_tilt_deg = np.rad2deg(refined_tilt_rad)

            # Rotate subtomograms in YX plane
            for _j, img in enumerate(subtomograms):
                img: ip.ImgArray
                angle = refined_tilt_deg[_j]
                img.rotate(-angle, cval=0, update=True)

            # zx-shift correction by self-ZNCC
            subtomo_proj = subtomograms.mean(axis="y")

            if edge_sigma is not None:
                # Regions outside the mask don't need to be considered.
                xc = int(subtomo_proj.shape.x / 2)
                w = int(spl.config.fit_width / scale / 2)
                subtomo_proj = subtomo_proj[ip.slicer.x[xc - w : xc + w + 1]]

            max_shift_px = max_shift / scale * 2
            pf_ang = 360 / spl.config.npf_range.center
            degrees = np.linspace(-pf_ang / 2, pf_ang / 2, n_rotations) + 180
            shifts = _misc.multi_rotated_auto_zncc(subtomo_proj, degrees, max_shift_px)

        # Update spline coordinates.
        # Because centers of subtomogram are on lattice points of pixel coordinate,
        # coordinates that will be shifted should be converted to integers.
        coords_px = self.nm2pixel(spl.map(anc), binsize=binsize).astype(np.float32)
        coords_px_new = _misc.shift_coords(coords_px, shifts, refined_tilt_rad)
        coords = coords_px_new * scale + self.multiscale_translation(binsize)

        # Update spline parameters
        self.splines[i] = spl.fit(coords, err_max=err_max)
        result = _misc.FitResult(shifts * scale)
        LOGGER.info(f" >> Shift RMSD = {result.rmsd:.3f} nm")
        return result

    @_misc.batch_process
    def fit_centroid(
        self,
        i: int = None,
        *,
        max_interval: nm = 30.0,
        binsize: int = 1,
        err_max: nm = 1.0,
        max_shift: nm = 5.0,
    ) -> _misc.FitResult:
        LOGGER.info(f"Running: {self.__class__.__name__}.fit_centroid, i={i}")
        spl = self.splines[i]
        anc = spl.prep_anchor_positions(max_interval=max_interval)
        scale = self.scale * binsize

        # sample subtomograms
        loader = _misc.prep_loader_for_refine(self, spl, anc, binsize)
        subtomograms = ip.asarray(loader.asnumpy(), axes="pzyx").mean(axis="y")[
            ip.slicer.x[::-1]
        ]
        num, lz, lx = subtomograms.shape
        dpx = ceilint(max_shift / scale)
        sl_z = slice(max((lz - 1) // 2 - dpx, 0), min(lz // 2 + dpx + 1, lz))
        sl_x = slice(max((lx - 1) // 2 - dpx, 0), min(lx // 2 + dpx + 1, lx))
        centers = np.stack(
            [centroid_2d(patch, sl_z, sl_x) for patch in subtomograms], axis=0
        )
        shifts = centers - np.column_stack(
            [
                np.full(num, (sl_z.start + sl_z.stop - 1) / 2),
                np.full(num, (sl_x.start + sl_x.stop - 1) / 2),
            ]
        )
        self.splines[i] = spl.shift(anc, shifts=shifts * scale, err_max=err_max)
        result = _misc.FitResult(shifts * scale)
        LOGGER.info(f" >> Shift RMSD = {result.rmsd:.3f} nm")
        return result

    @_misc.batch_process
    def refine(
        self,
        i: int = None,
        *,
        max_interval: nm = 30.0,
        binsize: int = 1,
        err_max: nm = 1.0,
        corr_allowed: float = 0.9,
        max_shift: nm = 2.0,
        n_rotations: int = 3,
    ) -> _misc.FitResult:
        """
        Spline refinement using global lattice structural parameters.

        Refine spline using the result of previous fit and the global structural parameters.
        During refinement, Y-projection of XZ cross section of cylinder is rotated with the
        twist angles, thus is much more precise than the coarse fitting.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to fit.
        max_interval : nm, default 24.0
            Maximum interval of sampling points in nm unit.
        binsize : int, default 1
            Multiscale bin size used for refining.
        corr_allowed : float, defaul is 0.9
            How many images will be used to make template for alignment. If 0.9, then top
            90% will be used.
        max_shift: nm, default 2.0
            Maximum shift from the true center of the cylinder. This parameter is used in
            phase cross correlation.
        n_rotations : int, default 3
            Number of rotations to be tested during finding the cylinder center.

        Returns
        -------
        FitResult
            Result of fitting.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.refine, i={i}")
        spl = self.splines[i]
        _required = [H.spacing, H.twist, H.npf]
        if not spl.props.has_glob(_required):
            if (radius := spl.radius) is None:
                radius = self.measure_radius(
                    i=i,
                    binsize=binsize,
                    positions="auto",
                    update=False,
                )
            with spl.props.temp_glob(radius=radius):
                gprops = self.global_cft_params(
                    i=i, binsize=binsize, nsamples=1, update=False
                )
        else:
            gprops = spl.props.glob.select(_required)
        gdict = {k: float(gprops[k][0]) for k in _required}
        ancs = spl.prep_anchor_positions(max_interval=max_interval)

        # Calculate Fourier parameters by cylindrical transformation along spline.
        # Skew angles are divided by the angle of single protofilament and the residual
        # angles are used, considering missing wedge effect.
        space = gdict[H.spacing]
        twist = gdict[H.twist]
        npf = roundint(gdict[H.npf])

        LOGGER.info(f" >> Parameters: spacing = {space:.2f} nm, twist = {twist:.3f} deg, PF = {npf}")  # fmt: skip

        # complement twisting
        pf_ang = 360 / npf
        twists = _misc.get_twists(spl.length(), ancs.size, space, twist, npf)
        scale = self.scale * binsize
        loader = _misc.prep_loader_for_refine(self, spl, ancs, binsize, twists)
        subtomograms = ip.asarray(loader.asnumpy(), axes="pzyx", dtype=np.float32)
        subtomograms[:] -= subtomograms.mean()  # normalize
        subtomograms.set_scale(zyx=scale)

        degrees = np.linspace(-pf_ang / 2, pf_ang / 2, n_rotations) + 180
        max_shift_px = max_shift / scale
        with set_gpu():
            inputs = subtomograms.mean(axis="y")[ip.slicer.x[::-1]]

            # Align twist-corrected images
            shifts_loc = _misc.multi_rotated_auto_zncc(inputs, degrees, max_shift_px)
            tasks = [
                _misc.delayed_translate(inputs[_j], shifts_loc[_j])
                for _j in range(ancs.size)
            ]
            imgs_aligned = _filter_by_corr(
                np.stack(compute(*tasks), axis="p"),
                corr_allowed,
            )

            # Make 2D template using coarse aligned images.
            imgcory = imgs_aligned.mean(axis="p")
            shift = rotated_auto_zncc(
                imgcory, degrees=degrees, max_shifts=max_shift_px * 2
            )
            template = imgcory.affine(translation=shift, mode=Mode.constant, cval=0.0)

            # Align twist-corrected images to the template
            quat = loader.molecules.quaternion()
            tasks = [
                _misc.delayed_zncc_maximum(
                    inputs[_j],
                    _misc.mask_missing_wedge(template, self.tilt_model, quat[_j]),
                    max_shift_px,
                    twists[_j],
                )
                for _j in range(ancs.size)
            ]
            shifts = np.stack(compute(*tasks), axis=0)

        # Update spline parameters
        self.splines[i] = spl.shift(ancs, shifts=shifts * scale, err_max=err_max)
        result = _misc.FitResult(shifts * scale)
        LOGGER.info(f" >> Shift RMSD = {result.rmsd:.3f} nm")
        return result

    @_misc.batch_process
    def measure_radius(
        self,
        i: int = None,
        *,
        binsize: int = 1,
        positions: NDArray[np.float32] | Literal["auto", "anchor"] = "auto",
        min_radius: nm = 1.0,
        max_radius: nm = 100.0,
        update: bool = True,
    ) -> nm:
        """
        Measure radius using radial profile from the center.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to measure.
        binsize : int, default 1
            Multiscale bin size used for radius calculation.
        positions : array-like or "auto" or "anchor", default "auto"
            Sampling positions (between 0 and 1) to calculate radius. If "anchor"
            is given, anchors of the spline will be used. If "auto" is given,
            three positions along the spline will be used.
        min_radius : nm, default 1.0
            Minimum radius of the cylinder.
        max_radius : nm, default 100.0
            Maximum radius of the cylinder.
        update : bool, default True
            If True, global radius property will be updated.

        Returns
        -------
        float (nm)
            Cylinder radius.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.measure_radius, i={i}")
        spl = self.splines[i]

        if isinstance(positions, str) and positions == "auto":
            nanchor = 3
            pos = 1 / nanchor * np.arange(nanchor) + 0.5 / nanchor
        elif isinstance(positions, str) and positions == "anchor":
            pos = spl.anchors
        else:
            pos = np.asarray(positions, dtype=np.float32)

        input_img = self._get_multiscale_or_original(binsize)

        depth = spl.config.fit_depth
        _scale = input_img.scale.x
        min_radius_px = min_radius / _scale
        max_radius = min(max_radius, spl.config.fit_width / 2)
        max_radius_px = max_radius / _scale
        spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
        tasks = [
            _misc.get_radial_prof(
                input_img, spl_trans, anc, (min_radius, max_radius), depth
            )
            for anc in pos
        ]
        profs: list[NDArray[np.float32]] = compute(*tasks)
        prof = np.stack(profs, axis=0).mean(axis=0)
        imax_sub = find_centroid_peak(prof, *_misc.get_thickness(spl, _scale))
        offset_px = _misc.get_radius_offset(min_radius_px, max_radius_px)
        radius = (imax_sub + offset_px) * _scale
        if update:
            spl.props.update_glob(
                [pl.Series(H.radius, [radius], dtype=pl.Float32)],
                bin_size=binsize,
            )
        LOGGER.info(f" >> Radius = {radius:.3f} nm")
        return radius

    @_misc.batch_process
    def local_radii(
        self,
        *,
        i: int = None,
        depth: nm = 50.0,
        binsize: int = 1,
        min_radius: nm = 1.0,
        max_radius: nm = 100.0,
        update: bool = True,
        update_glob: bool = True,
    ) -> pl.Series:
        """
        Measure the local radii along the splines.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        depth : nm, default 50.0
            Longitudinal length of subtomograms for calculation.
        binsize : int, default 1
            Multiscale binsize to be used.
        min_radius : nm, default 1.0
            Minimum radius of the cylinder.
        max_radius : nm, default 100.0
            Maximum radius of the cylinder.
        update : bool, default True
            If True, spline properties will be updated.
        update_glob : bool, default True
            If True, global properties will be updated using the mean of the local
            radii.

        Returns
        -------
        pl.Series
            Radii along the spline.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.local_radii, i={i}")
        spl = self.splines[i]

        input_img = self._get_multiscale_or_original(binsize)

        depth = spl.config.fit_depth
        _scale = input_img.scale.x
        thickness = _misc.get_thickness(spl, _scale)
        min_radius_px = min_radius / _scale
        max_radius = min(max_radius, spl.config.fit_width / 2)
        max_radius_px = max_radius / _scale
        offset_px = _misc.get_radius_offset(min_radius_px, max_radius_px)
        spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
        tasks = []
        for anc in spl_trans.anchors:
            task = _misc.get_radial_prof(
                input_img, spl_trans, anc, (min_radius, max_radius), depth
            )
            tasks.append(task)
        profs: list[NDArray[np.float32]] = compute(*tasks)
        radii = list[float]()
        for prof in profs:
            imax_sub = find_centroid_peak(prof, *thickness)
            radii.append((imax_sub + offset_px) * _scale)

        out = pl.Series(H.radius, radii, dtype=pl.Float32)
        if update:
            spl.props.update_loc([out], depth, bin_size=binsize)
        if update_glob:
            spl.props.update_glob([pl.Series(H.radius, [out.mean()])], bin_size=binsize)
        return out

    @_misc.batch_process
    def local_cft_params(
        self,
        *,
        i: int = None,
        depth: nm = 50.0,
        binsize: int = 1,
        radius: nm | Literal["local", "global"] = "global",
        nsamples: int = 8,
        update: bool = True,
        update_glob: bool = False,
    ) -> pl.DataFrame:
        """
        Calculate local lattice parameters from cylindrical Fourier space.

        To determine the peaks upsampled discrete Fourier transformation is used
        for every subtomogram.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        depth : nm, default 50.0
            Length of subtomogram for calculation of local parameters.
        binsize : int, default 1
            Multiscale bin size used for calculation.
        radius : str, default "global"
            If "local", use the local radius for the analysis. If "global", use the
            global radius. If a float, use the given radius.
        nsamples : int, default 8
            Number of cylindrical coordinate samplings for Fourier transformation. Multiple
            samplings are needed because up-sampled discrete Fourier transformation does not
            return exactly the same power spectra with shifted inputs, unlike FFT. Larger
            ``nsamples`` reduces the error but is slower.
        update : bool, default True
            If True, spline properties will be updated.
        update_glob : bool, default False
            If True, global properties will be updated using the mean or mode of the local
            properties.

        Returns
        -------
        polars.DataFrame
            Local properties.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.local_cft_params, i={i}")
        spl = self.splines[i]
        radii = _misc.prepare_radii(spl, radius)
        input_img = self._get_multiscale_or_original(binsize)
        _scale = input_img.scale.x
        tasks: list[Delayed[LatticeParams]] = []
        spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
        _analyze_fn = LatticeAnalyzer(spl.config).estimate_lattice_params_task
        for anc, r0 in zip(spl_trans.anchors, radii, strict=True):
            rmin, rmax = spl.radius_range(r0)
            rc = (rmin + rmax) / 2
            coords = spl_trans.local_cylindrical((rmin, rmax), depth, anc, scale=_scale)
            tasks.append(_analyze_fn(input_img, coords, rc, nsamples=nsamples))

        lprops = pl.DataFrame(
            compute(*tasks),
            schema=LatticeParams.polars_schema(),
        )
        if update:
            spl.props.update_loc(lprops, depth, bin_size=binsize)
        if update_glob:
            gprops = lprops.select(
                pl.col(H.spacing).mean(),
                pl.col(H.pitch).mean(),
                pl.col(H.twist).mean(),
                pl.col(H.skew).mean(),
                pl.col(H.rise).mean(),
                pl.col(H.rise_length).mean(),
                pl.col(H.npf).mode().first(),
                pl.col(H.start).mode().first(),
            )
            spl.props.update_glob(gprops, bin_size=binsize)

        return lprops

    def iter_local_image(
        self,
        i: int,
        depth: nm = 50.0,
        pos: int | None = None,
        binsize: int = 1,
    ) -> Iterable[ip.ImgArray]:
        spl = self.splines[i]
        if spl.radius is None:
            raise ValueError("Radius has not been determined yet.")

        input_img = self._get_multiscale_or_original(binsize)
        _scale = input_img.scale.x
        rmin, rmax = spl.radius_range()
        rc = (rmin + rmax) / 2
        if pos is None:
            anchors = spl.anchors
        else:
            anchors = [spl.anchors[pos]]
        spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
        for anc in anchors:
            coords = spl_trans.local_cylindrical((rmin, rmax), depth, anc, scale=_scale)
            polar = get_polar_image(input_img, coords, rc)
            polar[:] -= np.mean(polar)
            yield polar

    @_misc.batch_process
    def local_cft(
        self,
        *,
        i: int = None,
        depth: nm = 50.0,
        pos: int | None = None,
        binsize: int = 1,
    ) -> ip.ImgArray:
        """
        Calculate non-upsampled local cylindric Fourier transormation along spline.

        This function does not up-sample.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        depth : nm, default 50.0
            Length of subtomogram for calculation of local parameters.
        pos : int, optional
            Only calculate at ``pos``-th anchor if given.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            FT images stacked along "p" axis.
        """
        out = list[ip.ImgArray]()
        with set_gpu():
            for polar in self.iter_local_image(i, depth, pos, binsize):
                out.append(polar.fft(dims="rya"))
        return np.stack(out, axis="p")

    @_misc.batch_process
    def local_cps(
        self,
        *,
        i: int = None,
        depth: nm = 50.0,
        pos: int | None = None,
        binsize: int = 1,
    ) -> ip.ImgArray:
        """
        Calculate non-upsampled local cylindric power spectra along spline.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        depth : nm, default 50.0
            Length of subtomogram for calculation of local parameters.
        pos : int, optional
            Only calculate at ``pos``-th anchor if given.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            FT images stacked along "p" axis.
        """
        cft = self.local_cft(i=i, depth=depth, pos=pos, binsize=binsize)
        return cft.real**2 + cft.imag**2

    @_misc.batch_process
    def local_image_with_peaks(
        self,
        *,
        i: int = None,
        binsize: int | None = None,
    ) -> list[_misc.ImageWithPeak]:
        spl = self.splines[i]
        depth = spl.props.window_size[H.twist]
        if binsize is None:
            binsize = spl.props.binsize_loc[H.twist]
        df_loc = spl.props.loc
        out = list[_misc.ImageWithPeak]()
        for j, polar_img in enumerate(self.iter_local_image(i, depth, binsize=binsize)):
            cparams = spl.cylinder_params(
                spacing=_misc.get_component(df_loc, H.spacing, j),
                pitch=_misc.get_component(df_loc, H.pitch, j),
                twist=_misc.get_component(df_loc, H.twist, j),
                skew=_misc.get_component(df_loc, H.skew, j),
                rise_angle=_misc.get_component(df_loc, H.rise, j),
                start=_misc.get_component(df_loc, H.start, j),
                npf=_misc.get_component(df_loc, H.npf, j),
            )
            analyzer = LatticeAnalyzer(spl.config)
            peakv, peakh = analyzer.params_to_peaks(polar_img[0], cparams)
            peakv = peakv.shift_to_center()
            peakh = peakh.shift_to_center()
            out.append(_misc.ImageWithPeak(polar_img, [peakv, peakh]))
        return out

    @_misc.batch_process
    def global_cps(
        self,
        *,
        i: int = None,
        binsize: int = 1,
    ) -> ip.ImgArray:
        """
        Calculate global cylindrical power spectra.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            Complex image.
        """
        cft = self.global_cft(i=i, binsize=binsize)
        return cft.real**2 + cft.imag**2

    @_misc.batch_process
    def global_image_with_peaks(
        self,
        *,
        i: int = None,
        binsize: int | None = None,
    ) -> _misc.ImageWithPeak:
        spl = self.splines[i]
        if binsize is None:
            binsize = spl.props.binsize_glob[H.twist]
        img_st = self.straighten_cylindric(i, binsize=binsize)
        img_st -= np.mean(img_st)
        cparams = spl.cylinder_params(
            spacing=spl.props.get_glob(H.spacing, default=None),
            pitch=spl.props.get_glob(H.pitch, default=None),
            twist=spl.props.get_glob(H.twist, default=None),
            skew=spl.props.get_glob(H.skew, default=None),
            rise_angle=spl.props.get_glob(H.rise, default=None),
            start=spl.props.get_glob(H.start, default=None),
            npf=spl.props.get_glob(H.npf, default=None),
        )
        analyzer = LatticeAnalyzer(spl.config)
        peakv, peakh = analyzer.params_to_peaks(img_st[0], cparams)
        peakv = peakv.shift_to_center()
        peakh = peakh.shift_to_center()
        return _misc.ImageWithPeak(img_st, [peakv, peakh])

    @_misc.batch_process
    def global_cft_params(
        self,
        *,
        i: int = None,
        binsize: int = 1,
        nsamples: int = 8,
        update: bool = True,
    ) -> pl.DataFrame:
        """
        Calculate global lattice parameters.

        This function transforms tomogram using cylindrical coordinate system along
        spline. This function calls ``straighten`` beforehand, so that Fourier space is
        distorted if the cylindrical structure is curved.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        binsize : int, default 1
            Multiscale bin size used for calculation.
        nsamples : int, default 8
            Number of cylindrical coordinate samplings for Fourier transformation.
            Multiple samplings are needed because up-sampled discrete Fourier
            transformation does not return exactly the same power spectra with shifted
            inputs, unlike FFT. Larger ``nsamples`` reduces the error but is slower.
        update : bool, default True
            If True, spline properties will be updated.

        Returns
        -------
        pl.DataFrame
            Global properties.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.global_cft_params, i={i}")
        spl = self.splines[i]
        rmin, rmax = spl.radius_range()
        img_st = self.straighten_cylindric(i, radii=(rmin, rmax), binsize=binsize)
        rc = (rmin + rmax) / 2
        analyzer = LatticeAnalyzer(spl.config)
        lparams = analyzer.estimate_lattice_params_polar(img_st, rc, nsamples=nsamples)
        out = lparams.to_polars()
        if update:
            spl.props.update_glob(spl.props.glob.with_columns(out), bin_size=binsize)
        return out

    @_misc.batch_process
    def global_cft(self, i: int = None, binsize: int = 1) -> ip.ImgArray:
        """
        Calculate global cylindrical fast Fourier tranformation.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            Complex image.
        """
        img_st = self.straighten_cylindric(i, binsize=binsize)
        img_st -= np.mean(img_st)
        return img_st.fft(dims="rya")

    @_misc.batch_process
    def infer_polarity(
        self,
        i: int = None,
        *,
        binsize: int = 1,
        depth: nm = 40,
        mask_freq: bool = True,
        update: bool = True,
    ) -> Ori:
        """
        Infer spline polarities using polar 2D image.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to analyze.
        binsize : int, default 1
            Multiscale bin size used for calculation.
        depth : nm, default 40.0
            Depth of images used to infer polarities.

        Returns
        -------
        Ori
            Orientation of corresponding splines.
        """
        LOGGER.info(f"Running: {self.__class__.__name__}.infer_polarity, i={i}")
        current_scale = self.scale * binsize

        if binsize > 1:
            imgb = self.get_multiscale(binsize)
        else:
            try:
                imgb = self.get_multiscale(1)
            except ValueError:
                imgb = self.image

        spl = self.splines[i]
        cfg = spl.config
        ori_clockwise = Ori(cfg.clockwise)
        ori_counterclockwise = Ori.invert(ori_clockwise, allow_none=False)
        if spl.radius is None:
            r_range = 0.5 * current_scale, cfg.fit_width / 2
        else:
            r_range = spl.radius_range()
        point = 0.5  # the sampling point
        coords = spl.local_cylindrical(r_range, depth, point, scale=current_scale)
        polar = get_polar_image(imgb, coords, spl.radius, order=1)
        if mask_freq:
            polar = LatticeAnalyzer(cfg).mask_spectra(polar)
        img_flat = polar.mean(axis="y")

        if (npf := spl.props.get_glob(H.npf, None)) is None:
            # if the global properties are already calculated, use it
            # otherwise, calculate the number of PFs from the power spectrum
            ft = img_flat.fft(shift=False, dims="ra")
            pw = ft.real**2 + ft.imag**2
            img_pw = np.mean(pw, axis=0)
            npf = np.argmax(img_pw[cfg.npf_range.asslice()]) + cfg.npf_range.min

        pw_peak = img_flat.local_power_spectra(
            key=ip.slicer.a[npf - 1 : npf + 2],
            upsample_factor=20,
            dims="ra",
        ).mean(axis="a")
        r_argmax = np.argmax(pw_peak)
        clkwise = r_argmax - (pw_peak.size + 1) // 2 <= 0
        ori = ori_clockwise if clkwise else ori_counterclockwise

        # logging
        _val = pw_peak[r_argmax]
        pw_non_peak = np.delete(pw_peak, r_argmax)
        _ave, _std = np.mean(pw_non_peak), np.std(pw_non_peak, ddof=1)
        LOGGER.info(f" >> polarity = {ori.name} (peak intensity={_val:.2g} compared to {_ave:.2g} ± {_std:.2g})")  # fmt: skip
        if update:
            spl.orientation = ori
        return ori

    @_misc.batch_process
    def straighten(
        self,
        i: int = None,
        *,
        size: nm | tuple[nm, nm] = None,
        range_: tuple[float, float] = (0.0, 1.0),
        chunk_length: nm | None = None,
        binsize: int = 1,
    ) -> ip.ImgArray:
        """
        Straightening by building curved coordinate system around splines. Currently
        Cartesian coordinate system and cylindrical coordinate system are supported.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to straighten.
        size : float (nm), optional
            Vertical/horizontal box size.
        range_ : tuple[float, float], default (0.0, 1.0)
            Range of spline domain.
        chunk_length : nm, optional
            If spline is longer than this, it will be first split into chunks,
            straightened respectively and all the straightened images are concatenated
            afterward, to avoid loading entire image into memory.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            Straightened image. If Cartesian coordinate system is used, it will have "zyx".
        """
        spl = self.splines[i]
        input_img = self._get_multiscale_or_original(binsize)
        chunk_length = _misc.normalize_chunk_length(input_img, chunk_length)
        return _straighten.straighten(input_img, spl, range_, size)

    @_misc.batch_process
    def straighten_cylindric(
        self,
        i: int = None,
        *,
        radii: tuple[nm, nm] | None = None,
        range_: tuple[float, float] = (0.0, 1.0),
        chunk_length: nm | None = None,
        binsize: int = 1,
    ) -> ip.ImgArray:
        """
        Straightening by building curved coordinate system around splines. Currently
        Cartesian coordinate system and cylindrical coordinate system are supported.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that you want to straighten.
        radii : tuple of float (nm), optional
            Lower/upper limit of radius.
        range_ : tuple[float, float], default (0.0, 1.0)
            Range of spline domain.
        chunk_length : nm, optional
            If spline is longer than this, it will be first split into chunks,
            straightened respectively and all the straightened images are concatenated
            afterward, to avoid loading entire image into memory.
        binsize : int, default 1
            Multiscale bin size used for calculation.

        Returns
        -------
        ip.ImgArray
            Straightened image. If Cartesian coordinate system is used, it will have "zyx".
        """
        spl = self.splines[i]
        input_img = self._get_multiscale_or_original(binsize)
        chunk_length = _misc.normalize_chunk_length(input_img, chunk_length)
        return _straighten.straighten_cylindric(input_img, spl, range_, radii)

    @_misc.batch_process
    def map_centers(
        self,
        i: int = None,
        *,
        interval: nm = 1.0,
        orientation: Ori | str | None = None,
        rotate_molecules: bool = True,
    ) -> Molecules:
        """
        Mapping molecules along the center of a cylinder.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that mapping will be calculated.
        interval : float (nm), optional
            Interval of molecules.
        rotate_molecules : bool, default True
            If True, twist the molecule orientations according to the spline twist.

        Returns
        -------
        Molecules
            Molecules object with mapped coordinates and angles.
        """
        spl = self.splines[i]
        u = spl.prep_anchor_positions(interval=interval)
        if rotate_molecules:
            spacing = spl.props.get_glob(H.spacing)
            twist = spl.props.get_glob(H.twist) / 2
            rotation = np.deg2rad(spl.distances(u) / spacing * twist)
        else:
            rotation = None
        mole = spl.anchors_to_molecules(u, rotation=rotation)
        if spl._need_rotation(orientation):
            mole = mole.rotate_by_rotvec_internal([np.pi, 0, 0])
        return mole

    def get_cylinder_model(
        self,
        i: int,
        offsets: tuple[float, float] = (0.0, 0.0),
        **kwargs,
    ) -> CylinderModel:  # fmt: skip
        """
        Return the cylinder model at the given spline ID.

        Parameters
        ----------
        i : int
            Spline ID from which model will be created.
        offsets : tuple of float, optional
            Offset of the model. See :meth:`map_monomers` for details.

        Returns
        -------
        CylinderModel
            The cylinder model.
        """
        return self.splines[i].cylinder_model(offsets=offsets, **kwargs)

    @_misc.batch_process
    def map_monomers(
        self,
        i: int = None,
        *,
        offsets: tuple[nm, float] | None = None,
        orientation: Ori | str | None = None,
        extensions: tuple[int, int] = (0, 0),
        **kwargs,
    ) -> Molecules:
        """
        Map monomers in a regular cylinder shape.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that mapping will be calculated.
        offsets : tuple of float, optional
            The offset of origin of oblique coordinate system to map monomers.
        orientation : Ori or str, optional
            Orientation of the y-axis of each molecule.

        Returns
        -------
        Molecules
            Object that represents monomer positions and angles.
        """
        model = self.get_cylinder_model(i, offsets=offsets, **kwargs)
        ny, na = model.shape
        ext0, ext1 = extensions
        if ny + ext0 + ext1 < 0:
            raise ValueError("The number of monomers is negative.")
        yy, aa = np.indices((ny + ext0 + ext1, na), dtype=np.int32)
        yy -= ext0
        coords = np.stack([yy.ravel(), aa.ravel()], axis=1)
        spl = self.splines[i]
        mole = model.locate_molecules(spl, coords)
        if spl._need_rotation(orientation):
            mole = mole.rotate_by_rotvec_internal([np.pi, 0, 0])
        return mole

    @_misc.batch_process
    def map_on_grid(
        self,
        i: int = None,
        coords: NDArray[np.int32] = (),
        *,
        offsets: tuple[nm, float] | None = None,
        orientation: Ori | str | None = None,
        **kwargs,
    ) -> Molecules:
        """
        Map monomers in a regular cylinder shape.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that mapping will be calculated.
        coords : ndarray
            Integer coordinates on the cylinder surface.
        offsets : tuple of float, optional
            The offset of origin of oblique coordinate system to map monomers.
        orientation : Ori or str, optional
            Orientation of the y-axis of each molecule.

        Returns
        -------
        Molecules
            Object that represents monomer positions and angles.
        """
        model = self.get_cylinder_model(i, offsets=offsets, **kwargs)
        coords = np.asarray(coords, dtype=np.int32)
        spl = self.splines[i]
        mole = model.locate_molecules(spl, coords)
        if spl._need_rotation(orientation):
            mole = mole.rotate_by_rotvec_internal([np.pi, 0, 0])
        return mole

    @_misc.batch_process
    def map_pf_line(
        self,
        i: int = None,
        *,
        interval: nm = 1.0,
        offsets: tuple[nm, float] = (0.0, 0.0),
        orientation: Ori | str | None = None,
    ) -> Molecules:
        """
        Mapping molecules along a protofilament line.

        This method is useful when you want to visualize seam or protofilament, or
        assign molecules for subtomogram averaging of seam binding protein or doublet
        microtubule.

        Parameters
        ----------
        i : int or iterable of int, optional
            Spline ID that mapping will be calculated.
        offsets : (float, float), default (0.0, 0.0)
            Axial offset in nm and angular offset in degree.

        Returns
        -------
        Molecules
            Object that represents protofilament positions and angles.
        """
        spl = self.splines[i]
        spacing = spl.props.get_glob(H.spacing)
        twist = spl.props.get_glob(H.twist) / 2

        ny = roundint(spl.length() / interval)
        skew_rad = np.deg2rad(twist) * interval / spacing

        yoffset, aoffset = offsets
        rcoords = np.full(ny, spl.radius)
        ycoords = np.arange(ny) * interval + yoffset
        acoords = np.arange(ny) * skew_rad + np.deg2rad(aoffset)
        coords = np.stack([rcoords, ycoords, acoords], axis=1)
        mole = spl.cylindrical_to_molecules(coords)
        if spl._need_rotation(orientation):
            mole = mole.rotate_by_rotvec_internal([np.pi, 0, 0])
        return mole

    def _prep_fit_spline(self, spl: CylSpline, anc: NDArray[np.float32], binsize: int):
        npoints = anc.size
        interval = spl.length() / (npoints - 1)
        depth_px = self.nm2pixel(spl.config.fit_depth, binsize=binsize)
        width_px = self.nm2pixel(spl.config.fit_width, binsize=binsize)

        # If subtomogram region is rotated by 45 degree, its XY-width will be
        # (length + width) / sqrt(2)
        if binsize > 1:
            centers = spl.map(anc) - self.multiscale_translation(binsize)
        else:
            centers = spl.map(anc)
        center_px = self.nm2pixel(centers, binsize=binsize)
        size_px = (width_px,) + (roundint((width_px + depth_px) / 1.414),) * 2
        input_img = self._get_multiscale_or_original(binsize)

        subtomograms = crop_tomograms(input_img, center_px, size_px)
        subtomograms[:] -= subtomograms.mean()
        scale = self.scale * binsize
        return subtomograms, interval, scale
splines property

List of splines.

add_spline(coords, *, order=3, err_max=0.5, extrapolate=ExtrapolationMode.linear, config={})

Add spline path to tomogram.

Parameters:

Name Type Description Default
coords array - like

(N, 3) array of coordinates. A spline curve that fit it well is added.

 required
order int

Order of spline curve.

 3
extrapolate str

Extrapolation mode of the spline.

 linear
config SplineConfig or dict

Configuration for spline fitting.

 {}
Source code in cylindra/components/tomogram/_cyl_tomo.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def add_spline(
    self,
    coords: ArrayLike,
    *,
    order: int = 3,
    err_max: nm = 0.5,
    extrapolate: ExtrapolationMode | str = ExtrapolationMode.linear,
    config: SplineConfig | dict[str, Any] = {},
) -> None:
    """
    Add spline path to tomogram.

    Parameters
    ----------
    coords : array-like
        (N, 3) array of coordinates. A spline curve that fit it well is added.
    order : int, optional
        Order of spline curve.
    extrapolate : str, optional
        Extrapolation mode of the spline.
    config : SplineConfig or dict, optional
        Configuration for spline fitting.
    """
    _coords = np.asarray(coords)
    ncoords = _coords.shape[0]
    spl = CylSpline(
        order=order,
        config=config,
        extrapolate=extrapolate,
    ).fit(_coords, err_max=err_max)
    interval: nm = 30.0
    length = spl.length()

    n = int(length / interval) + 1
    fit = spl.map(np.linspace(0, 1, n))
    if ncoords <= spl.order and ncoords < fit.shape[0]:
        return self.add_spline(
            fit,
            order=order,
            err_max=err_max,
            extrapolate=extrapolate,
            config=config,
        )

    self.splines.append(spl)
    return None
align_to_polarity(orientation=Ori.MinusToPlus)

Align all the splines in the direction parallel to the given polarity.

Parameters:

Name Type Description Default
orientation Ori or str

To which direction splines will be aligned.

 Ori.MinusToPlus

Returns:

Type Description
Tomogram object

Same object with updated splines.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def align_to_polarity(self, orientation: Ori | str = Ori.MinusToPlus) -> Self:
    """
    Align all the splines in the direction parallel to the given polarity.

    Parameters
    ----------
    orientation : Ori or str, default Ori.MinusToPlus
        To which direction splines will be aligned.

    Returns
    -------
    Tomogram object
        Same object with updated splines.
    """
    orientation = Ori(orientation)
    if orientation is Ori.none:
        raise ValueError("Must be PlusToMinus or MinusToPlus.")
    for i, spl in enumerate(self.splines):
        if spl.orientation is Ori.none:
            raise ValueError(f"Spline-{i} has no orientation.")
        if spl.orientation != orientation:
            try:
                self.splines[i] = spl.invert()
            except Exception as e:
                raise type(e)(f"Cannot invert spline-{i}: {e}")
    return self
fit(i=None, *, max_interval=30.0, degree_precision=0.5, binsize=1, err_max=1.0, edge_sigma=2.0, max_shift=5.0, n_rotations=5)

Roughly fit splines to cylindrical structures.

Subtomograms will be sampled at every max_interval nm. In dense mode, Subtomograms will be masked relative to XY-plane, using sigmoid function. Sharpness of the sigmoid function is determined by dense_mode_sigma (dense_mode_sigma=0 corresponds to a step function).

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to fit.

 None
max_interval nm

Maximum interval of sampling points in nm unit.

 30.0
degree_precision float

Precision of xy-tilt degree in angular correlation.

 0.5
binsize int

Multiscale bin size used for fitting.

 1
edge_sigma nm

Sharpness of mask at the edges. If not None, fitting will be executed after regions outside the cylinder are masked. Soft mask is important for precision because sharp changes in intensity cause strong correlation at the edges.

 2.0
max_shift nm

Maximum shift from the true center of the cylinder. This parameter is used in phase cross correlation.

 5.0
n_rotations int

Number of rotations to be tested during finding the cylinder center.

 5

Returns:

Type Description
FitResult

Result of fitting.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
@_misc.batch_process
def fit(
    self,
    i: int = None,
    *,
    max_interval: nm = 30.0,
    degree_precision: float = 0.5,
    binsize: int = 1,
    err_max: nm = 1.0,
    edge_sigma: nm = 2.0,
    max_shift: nm = 5.0,
    n_rotations: int = 5,
) -> _misc.FitResult:
    """
    Roughly fit splines to cylindrical structures.

    Subtomograms will be sampled at every ``max_interval`` nm. In dense mode,
    Subtomograms will be masked relative to XY-plane, using sigmoid function.
    Sharpness of the sigmoid function is determined by ``dense_mode_sigma``
    (``dense_mode_sigma=0`` corresponds to a step function).

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to fit.
    max_interval : nm, default 30.0
        Maximum interval of sampling points in nm unit.
    degree_precision : float, default 0.5
        Precision of xy-tilt degree in angular correlation.
    binsize : int, default 1
        Multiscale bin size used for fitting.
    edge_sigma : nm, default 2.0
        Sharpness of mask at the edges. If not None, fitting will be executed after regions
        outside the cylinder are masked. Soft mask is important for precision because sharp
        changes in intensity cause strong correlation at the edges.
    max_shift: nm, default 5.0
        Maximum shift from the true center of the cylinder. This parameter is used in phase
        cross correlation.
    n_rotations : int, default 5
        Number of rotations to be tested during finding the cylinder center.

    Returns
    -------
    FitResult
        Result of fitting.
    """
    LOGGER.info(f"Running: {self.__class__.__name__}.fit, i={i}")
    spl = self.splines[i]
    anc = spl.prep_anchor_positions(max_interval=max_interval)
    subtomograms, interval, scale = self._prep_fit_spline(spl, anc, binsize)

    with set_gpu():
        subtomograms = _misc.soft_mask_edges(
            subtomograms, spl, anc, scale, edge_sigma
        )
        ds = spl.map(anc, der=1)
        yx_tilt = np.rad2deg(np.arctan2(-ds[:, 2], ds[:, 1]))
        degree_max = 14.0
        nrots = roundint(degree_max / degree_precision) + 1

        # Angular correlation
        out = _misc.dask_angle_corr(subtomograms, yx_tilt, nrots=nrots)
        refined_tilt_deg = np.array(out)
        refined_tilt_rad = np.deg2rad(refined_tilt_deg)

        # If subtomograms are sampled at short intervals, angles should be smoothened to
        # avoid overfitting.
        size = 2 * roundint(48.0 / interval) + 1
        if size > 1:
            # Mirror-mode padding is "a b c d | c b a".
            refined_tilt_rad = _misc.angle_uniform_filter(
                refined_tilt_rad, size=size, mode=Mode.mirror
            )
            refined_tilt_deg = np.rad2deg(refined_tilt_rad)

        # Rotate subtomograms in YX plane
        for _j, img in enumerate(subtomograms):
            img: ip.ImgArray
            angle = refined_tilt_deg[_j]
            img.rotate(-angle, cval=0, update=True)

        # zx-shift correction by self-ZNCC
        subtomo_proj = subtomograms.mean(axis="y")

        if edge_sigma is not None:
            # Regions outside the mask don't need to be considered.
            xc = int(subtomo_proj.shape.x / 2)
            w = int(spl.config.fit_width / scale / 2)
            subtomo_proj = subtomo_proj[ip.slicer.x[xc - w : xc + w + 1]]

        max_shift_px = max_shift / scale * 2
        pf_ang = 360 / spl.config.npf_range.center
        degrees = np.linspace(-pf_ang / 2, pf_ang / 2, n_rotations) + 180
        shifts = _misc.multi_rotated_auto_zncc(subtomo_proj, degrees, max_shift_px)

    # Update spline coordinates.
    # Because centers of subtomogram are on lattice points of pixel coordinate,
    # coordinates that will be shifted should be converted to integers.
    coords_px = self.nm2pixel(spl.map(anc), binsize=binsize).astype(np.float32)
    coords_px_new = _misc.shift_coords(coords_px, shifts, refined_tilt_rad)
    coords = coords_px_new * scale + self.multiscale_translation(binsize)

    # Update spline parameters
    self.splines[i] = spl.fit(coords, err_max=err_max)
    result = _misc.FitResult(shifts * scale)
    LOGGER.info(f" >> Shift RMSD = {result.rmsd:.3f} nm")
    return result
get_cylinder_model(i, offsets=(0.0, 0.0), **kwargs)

Return the cylinder model at the given spline ID.

Parameters:

Name Type Description Default
i int

Spline ID from which model will be created.

 required
offsets tuple of float

Offset of the model. See :meth:map_monomers for details.

 (0.0, 0.0)

Returns:

Type Description
CylinderModel

The cylinder model.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
def get_cylinder_model(
    self,
    i: int,
    offsets: tuple[float, float] = (0.0, 0.0),
    **kwargs,
) -> CylinderModel:  # fmt: skip
    """
    Return the cylinder model at the given spline ID.

    Parameters
    ----------
    i : int
        Spline ID from which model will be created.
    offsets : tuple of float, optional
        Offset of the model. See :meth:`map_monomers` for details.

    Returns
    -------
    CylinderModel
        The cylinder model.
    """
    return self.splines[i].cylinder_model(offsets=offsets, **kwargs)
global_cft(i=None, binsize=1)

Calculate global cylindrical fast Fourier tranformation.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
binsize int

Multiscale bin size used for calculation.

 1

Returns:

Type Description
ImgArray

Complex image.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
@_misc.batch_process
def global_cft(self, i: int = None, binsize: int = 1) -> ip.ImgArray:
    """
    Calculate global cylindrical fast Fourier tranformation.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    binsize : int, default 1
        Multiscale bin size used for calculation.

    Returns
    -------
    ip.ImgArray
        Complex image.
    """
    img_st = self.straighten_cylindric(i, binsize=binsize)
    img_st -= np.mean(img_st)
    return img_st.fft(dims="rya")
global_cft_params(*, i=None, binsize=1, nsamples=8, update=True)

Calculate global lattice parameters.

This function transforms tomogram using cylindrical coordinate system along spline. This function calls straighten beforehand, so that Fourier space is distorted if the cylindrical structure is curved.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
binsize int

Multiscale bin size used for calculation.

 1
nsamples int

Number of cylindrical coordinate samplings for Fourier transformation. Multiple samplings are needed because up-sampled discrete Fourier transformation does not return exactly the same power spectra with shifted inputs, unlike FFT. Larger nsamples reduces the error but is slower.

 8
update bool

If True, spline properties will be updated.

 True

Returns:

Type Description
DataFrame

Global properties.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
@_misc.batch_process
def global_cft_params(
    self,
    *,
    i: int = None,
    binsize: int = 1,
    nsamples: int = 8,
    update: bool = True,
) -> pl.DataFrame:
    """
    Calculate global lattice parameters.

    This function transforms tomogram using cylindrical coordinate system along
    spline. This function calls ``straighten`` beforehand, so that Fourier space is
    distorted if the cylindrical structure is curved.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    binsize : int, default 1
        Multiscale bin size used for calculation.
    nsamples : int, default 8
        Number of cylindrical coordinate samplings for Fourier transformation.
        Multiple samplings are needed because up-sampled discrete Fourier
        transformation does not return exactly the same power spectra with shifted
        inputs, unlike FFT. Larger ``nsamples`` reduces the error but is slower.
    update : bool, default True
        If True, spline properties will be updated.

    Returns
    -------
    pl.DataFrame
        Global properties.
    """
    LOGGER.info(f"Running: {self.__class__.__name__}.global_cft_params, i={i}")
    spl = self.splines[i]
    rmin, rmax = spl.radius_range()
    img_st = self.straighten_cylindric(i, radii=(rmin, rmax), binsize=binsize)
    rc = (rmin + rmax) / 2
    analyzer = LatticeAnalyzer(spl.config)
    lparams = analyzer.estimate_lattice_params_polar(img_st, rc, nsamples=nsamples)
    out = lparams.to_polars()
    if update:
        spl.props.update_glob(spl.props.glob.with_columns(out), bin_size=binsize)
    return out
global_cps(*, i=None, binsize=1)

Calculate global cylindrical power spectra.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
binsize int

Multiscale bin size used for calculation.

 1

Returns:

Type Description
ImgArray

Complex image.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
@_misc.batch_process
def global_cps(
    self,
    *,
    i: int = None,
    binsize: int = 1,
) -> ip.ImgArray:
    """
    Calculate global cylindrical power spectra.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    binsize : int, default 1
        Multiscale bin size used for calculation.

    Returns
    -------
    ip.ImgArray
        Complex image.
    """
    cft = self.global_cft(i=i, binsize=binsize)
    return cft.real**2 + cft.imag**2
infer_polarity(i=None, *, binsize=1, depth=40, mask_freq=True, update=True)

Infer spline polarities using polar 2D image.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
binsize int

Multiscale bin size used for calculation.

 1
depth nm

Depth of images used to infer polarities.

 40.0

Returns:

Type Description
Ori

Orientation of corresponding splines.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
@_misc.batch_process
def infer_polarity(
    self,
    i: int = None,
    *,
    binsize: int = 1,
    depth: nm = 40,
    mask_freq: bool = True,
    update: bool = True,
) -> Ori:
    """
    Infer spline polarities using polar 2D image.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    binsize : int, default 1
        Multiscale bin size used for calculation.
    depth : nm, default 40.0
        Depth of images used to infer polarities.

    Returns
    -------
    Ori
        Orientation of corresponding splines.
    """
    LOGGER.info(f"Running: {self.__class__.__name__}.infer_polarity, i={i}")
    current_scale = self.scale * binsize

    if binsize > 1:
        imgb = self.get_multiscale(binsize)
    else:
        try:
            imgb = self.get_multiscale(1)
        except ValueError:
            imgb = self.image

    spl = self.splines[i]
    cfg = spl.config
    ori_clockwise = Ori(cfg.clockwise)
    ori_counterclockwise = Ori.invert(ori_clockwise, allow_none=False)
    if spl.radius is None:
        r_range = 0.5 * current_scale, cfg.fit_width / 2
    else:
        r_range = spl.radius_range()
    point = 0.5  # the sampling point
    coords = spl.local_cylindrical(r_range, depth, point, scale=current_scale)
    polar = get_polar_image(imgb, coords, spl.radius, order=1)
    if mask_freq:
        polar = LatticeAnalyzer(cfg).mask_spectra(polar)
    img_flat = polar.mean(axis="y")

    if (npf := spl.props.get_glob(H.npf, None)) is None:
        # if the global properties are already calculated, use it
        # otherwise, calculate the number of PFs from the power spectrum
        ft = img_flat.fft(shift=False, dims="ra")
        pw = ft.real**2 + ft.imag**2
        img_pw = np.mean(pw, axis=0)
        npf = np.argmax(img_pw[cfg.npf_range.asslice()]) + cfg.npf_range.min

    pw_peak = img_flat.local_power_spectra(
        key=ip.slicer.a[npf - 1 : npf + 2],
        upsample_factor=20,
        dims="ra",
    ).mean(axis="a")
    r_argmax = np.argmax(pw_peak)
    clkwise = r_argmax - (pw_peak.size + 1) // 2 <= 0
    ori = ori_clockwise if clkwise else ori_counterclockwise

    # logging
    _val = pw_peak[r_argmax]
    pw_non_peak = np.delete(pw_peak, r_argmax)
    _ave, _std = np.mean(pw_non_peak), np.std(pw_non_peak, ddof=1)
    LOGGER.info(f" >> polarity = {ori.name} (peak intensity={_val:.2g} compared to {_ave:.2g} ± {_std:.2g})")  # fmt: skip
    if update:
        spl.orientation = ori
    return ori
local_cft(*, i=None, depth=50.0, pos=None, binsize=1)

Calculate non-upsampled local cylindric Fourier transormation along spline.

This function does not up-sample.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
depth nm

Length of subtomogram for calculation of local parameters.

 50.0
pos int

Only calculate at pos-th anchor if given.

 None
binsize int

Multiscale bin size used for calculation.

 1

Returns:

Type Description
ImgArray

FT images stacked along "p" axis.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
@_misc.batch_process
def local_cft(
    self,
    *,
    i: int = None,
    depth: nm = 50.0,
    pos: int | None = None,
    binsize: int = 1,
) -> ip.ImgArray:
    """
    Calculate non-upsampled local cylindric Fourier transormation along spline.

    This function does not up-sample.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    depth : nm, default 50.0
        Length of subtomogram for calculation of local parameters.
    pos : int, optional
        Only calculate at ``pos``-th anchor if given.
    binsize : int, default 1
        Multiscale bin size used for calculation.

    Returns
    -------
    ip.ImgArray
        FT images stacked along "p" axis.
    """
    out = list[ip.ImgArray]()
    with set_gpu():
        for polar in self.iter_local_image(i, depth, pos, binsize):
            out.append(polar.fft(dims="rya"))
    return np.stack(out, axis="p")
local_cft_params(*, i=None, depth=50.0, binsize=1, radius='global', nsamples=8, update=True, update_glob=False)

Calculate local lattice parameters from cylindrical Fourier space.

To determine the peaks upsampled discrete Fourier transformation is used for every subtomogram.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
depth nm

Length of subtomogram for calculation of local parameters.

 50.0
binsize int

Multiscale bin size used for calculation.

 1
radius str

If "local", use the local radius for the analysis. If "global", use the global radius. If a float, use the given radius.

 "global"
nsamples int

Number of cylindrical coordinate samplings for Fourier transformation. Multiple samplings are needed because up-sampled discrete Fourier transformation does not return exactly the same power spectra with shifted inputs, unlike FFT. Larger nsamples reduces the error but is slower.

 8
update bool

If True, spline properties will be updated.

 True
update_glob bool

If True, global properties will be updated using the mean or mode of the local properties.

 False

Returns:

Type Description
DataFrame

Local properties.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
@_misc.batch_process
def local_cft_params(
    self,
    *,
    i: int = None,
    depth: nm = 50.0,
    binsize: int = 1,
    radius: nm | Literal["local", "global"] = "global",
    nsamples: int = 8,
    update: bool = True,
    update_glob: bool = False,
) -> pl.DataFrame:
    """
    Calculate local lattice parameters from cylindrical Fourier space.

    To determine the peaks upsampled discrete Fourier transformation is used
    for every subtomogram.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    depth : nm, default 50.0
        Length of subtomogram for calculation of local parameters.
    binsize : int, default 1
        Multiscale bin size used for calculation.
    radius : str, default "global"
        If "local", use the local radius for the analysis. If "global", use the
        global radius. If a float, use the given radius.
    nsamples : int, default 8
        Number of cylindrical coordinate samplings for Fourier transformation. Multiple
        samplings are needed because up-sampled discrete Fourier transformation does not
        return exactly the same power spectra with shifted inputs, unlike FFT. Larger
        ``nsamples`` reduces the error but is slower.
    update : bool, default True
        If True, spline properties will be updated.
    update_glob : bool, default False
        If True, global properties will be updated using the mean or mode of the local
        properties.

    Returns
    -------
    polars.DataFrame
        Local properties.
    """
    LOGGER.info(f"Running: {self.__class__.__name__}.local_cft_params, i={i}")
    spl = self.splines[i]
    radii = _misc.prepare_radii(spl, radius)
    input_img = self._get_multiscale_or_original(binsize)
    _scale = input_img.scale.x
    tasks: list[Delayed[LatticeParams]] = []
    spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
    _analyze_fn = LatticeAnalyzer(spl.config).estimate_lattice_params_task
    for anc, r0 in zip(spl_trans.anchors, radii, strict=True):
        rmin, rmax = spl.radius_range(r0)
        rc = (rmin + rmax) / 2
        coords = spl_trans.local_cylindrical((rmin, rmax), depth, anc, scale=_scale)
        tasks.append(_analyze_fn(input_img, coords, rc, nsamples=nsamples))

    lprops = pl.DataFrame(
        compute(*tasks),
        schema=LatticeParams.polars_schema(),
    )
    if update:
        spl.props.update_loc(lprops, depth, bin_size=binsize)
    if update_glob:
        gprops = lprops.select(
            pl.col(H.spacing).mean(),
            pl.col(H.pitch).mean(),
            pl.col(H.twist).mean(),
            pl.col(H.skew).mean(),
            pl.col(H.rise).mean(),
            pl.col(H.rise_length).mean(),
            pl.col(H.npf).mode().first(),
            pl.col(H.start).mode().first(),
        )
        spl.props.update_glob(gprops, bin_size=binsize)

    return lprops
local_cps(*, i=None, depth=50.0, pos=None, binsize=1)

Calculate non-upsampled local cylindric power spectra along spline.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
depth nm

Length of subtomogram for calculation of local parameters.

 50.0
pos int

Only calculate at pos-th anchor if given.

 None
binsize int

Multiscale bin size used for calculation.

 1

Returns:

Type Description
ImgArray

FT images stacked along "p" axis.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
@_misc.batch_process
def local_cps(
    self,
    *,
    i: int = None,
    depth: nm = 50.0,
    pos: int | None = None,
    binsize: int = 1,
) -> ip.ImgArray:
    """
    Calculate non-upsampled local cylindric power spectra along spline.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    depth : nm, default 50.0
        Length of subtomogram for calculation of local parameters.
    pos : int, optional
        Only calculate at ``pos``-th anchor if given.
    binsize : int, default 1
        Multiscale bin size used for calculation.

    Returns
    -------
    ip.ImgArray
        FT images stacked along "p" axis.
    """
    cft = self.local_cft(i=i, depth=depth, pos=pos, binsize=binsize)
    return cft.real**2 + cft.imag**2
local_radii(*, i=None, depth=50.0, binsize=1, min_radius=1.0, max_radius=100.0, update=True, update_glob=True)

Measure the local radii along the splines.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that you want to analyze.

 None
depth nm

Longitudinal length of subtomograms for calculation.

 50.0
binsize int

Multiscale binsize to be used.

 1
min_radius nm

Minimum radius of the cylinder.

 1.0
max_radius nm

Maximum radius of the cylinder.

 100.0
update bool

If True, spline properties will be updated.

 True
update_glob bool

If True, global properties will be updated using the mean of the local radii.

 True

Returns:

Type Description
Series

Radii along the spline.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
@_misc.batch_process
def local_radii(
    self,
    *,
    i: int = None,
    depth: nm = 50.0,
    binsize: int = 1,
    min_radius: nm = 1.0,
    max_radius: nm = 100.0,
    update: bool = True,
    update_glob: bool = True,
) -> pl.Series:
    """
    Measure the local radii along the splines.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that you want to analyze.
    depth : nm, default 50.0
        Longitudinal length of subtomograms for calculation.
    binsize : int, default 1
        Multiscale binsize to be used.
    min_radius : nm, default 1.0
        Minimum radius of the cylinder.
    max_radius : nm, default 100.0
        Maximum radius of the cylinder.
    update : bool, default True
        If True, spline properties will be updated.
    update_glob : bool, default True
        If True, global properties will be updated using the mean of the local
        radii.

    Returns
    -------
    pl.Series
        Radii along the spline.
    """
    LOGGER.info(f"Running: {self.__class__.__name__}.local_radii, i={i}")
    spl = self.splines[i]

    input_img = self._get_multiscale_or_original(binsize)

    depth = spl.config.fit_depth
    _scale = input_img.scale.x
    thickness = _misc.get_thickness(spl, _scale)
    min_radius_px = min_radius / _scale
    max_radius = min(max_radius, spl.config.fit_width / 2)
    max_radius_px = max_radius / _scale
    offset_px = _misc.get_radius_offset(min_radius_px, max_radius_px)
    spl_trans = spl.translate([-self.multiscale_translation(binsize)] * 3)
    tasks = []
    for anc in spl_trans.anchors:
        task = _misc.get_radial_prof(
            input_img, spl_trans, anc, (min_radius, max_radius), depth
        )
        tasks.append(task)
    profs: list[NDArray[np.float32]] = compute(*tasks)
    radii = list[float]()
    for prof in profs:
        imax_sub = find_centroid_peak(prof, *thickness)
        radii.append((imax_sub + offset_px) * _scale)

    out = pl.Series(H.radius, radii, dtype=pl.Float32)
    if update:
        spl.props.update_loc([out], depth, bin_size=binsize)
    if update_glob:
        spl.props.update_glob([pl.Series(H.radius, [out.mean()])], bin_size=binsize)
    return out
make_anchors(i=None, *, interval=None, n=None, max_interval=None)

Make anchors on spline object(s).

Parameters:

Name Type Description Default
interval nm

Anchor intervals.

 None
n int

Number of anchors

 None
max_interval nm

Maximum interval between anchors.

 None
Source code in cylindra/components/tomogram/_cyl_tomo.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
@_misc.batch_process
def make_anchors(
    self,
    i: int = None,
    *,
    interval: nm | None = None,
    n: int | None = None,
    max_interval: nm | None = None,
):
    """
    Make anchors on spline object(s).

    Parameters
    ----------
    interval : nm, optional
        Anchor intervals.
    n : int, optional
        Number of anchors
    max_interval : nm, optional
        Maximum interval between anchors.

    """
    self.splines[i].make_anchors(interval=interval, n=n, max_interval=max_interval)
    return None
map_centers(i=None, *, interval=1.0, orientation=None, rotate_molecules=True)

Mapping molecules along the center of a cylinder.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that mapping will be calculated.

 None
interval float(nm)

Interval of molecules.

 1.0
rotate_molecules bool

If True, twist the molecule orientations according to the spline twist.

 True

Returns:

Type Description
Molecules

Molecules object with mapped coordinates and angles.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
@_misc.batch_process
def map_centers(
    self,
    i: int = None,
    *,
    interval: nm = 1.0,
    orientation: Ori | str | None = None,
    rotate_molecules: bool = True,
) -> Molecules:
    """
    Mapping molecules along the center of a cylinder.

    Parameters
    ----------
    i : int or iterable of int, optional
        Spline ID that mapping will be calculated.
    interval : float (nm), optional
        Interval of molecules.
    rotate_molecules : bool, default True
        If True, twist the molecule orientations according to the spline twist.

    Returns
    -------
    Molecules
        Molecules object with mapped coordinates and angles.
    """
    spl = self.splines[i]
    u = spl.prep_anchor_positions(interval=interval)
    if rotate_molecules:
        spacing = spl.props.get_glob(H.spacing)
        twist = spl.props.get_glob(H.twist) / 2
        rotation = np.deg2rad(spl.distances(u) / spacing * twist)
    else:
        rotation = None
    mole = spl.anchors_to_molecules(u, rotation=rotation)
    if spl._need_rotation(orientation):
        mole = mole.rotate_by_rotvec_internal([np.pi, 0, 0])
    return mole
map_monomers(i=None, *, offsets=None, orientation=None, extensions=(0, 0), **kwargs)

Map monomers in a regular cylinder shape.

Parameters:

Name Type Description Default
i int or iterable of int

Spline ID that mapping will be calculated.

 None
offsets tuple of float

The offset of origin of oblique coordinate system to map monomers.

 None
orientation Ori or str

Orientation of the y-axis of each molecule.

 None

Returns:

Type Description
Molecules

Object that represents monomer positions and angles.
Source code in cylindra/components/tomogram/_cyl_tomo.py 
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120