susan.utils.functions

susan.utils._functions.dose_from_fsc(fsc, apix, freq_range=(0.1, 0.8), fsc_min=0.1)

Estimate effective dose from the Guinier slope of the FSC curve.

The ExpFilt dose is applied in reconstruction as exp(−s²·dose/4), where s is in 1/Å. In the intermediate frequency range the FSC decays as the same Gaussian envelope, so fitting ln(FSC) vs s² gives slope = −dose/4, and:

dose = −4 · d(ln FSC)/d(s²)

This can be compared to the mean of ptcls.def_ExFl (excluding failures marked as 9999) to calibrate aligner.expfilt_gain:

expfilt_gain = dose_from_fsc(fsc, apix) / mean_estimated_dose

Parameters:

• fsc (array_like) – FSC curve as returned by fsc_get. Assumed to have n shells spanning a box of side 2(n-1), i.e. shell k → s = k / (2(n-1)·apix), matching fsc_get and fsc_analyse.

• apix (float) – Pixel size in Angstroms.

• freq_range (tuple of float) – (low, high) as fractions of Nyquist over which to fit. The default (0.1, 0.8) covers the Guinier decay while stopping before the noise-dominated tail.

• fsc_min (float) – Minimum FSC value included in the fit. Shells at or below the noise floor would bias the slope. Default 0.1.

Returns:

Effective dose in Ų consistent with the ExpFilt convention. Returns NaN if the fit cannot be performed.

Return type:

float

susan.utils._functions.radial_average(v)

Compute the radial (shell) average of a 3-D volume.

Each output bin k contains the mean of all voxels at radius r ≈ k pixels from the volume centre. The output length is set by the largest dimension so that all voxels are included.

Parameters:

v (ndarray, shape (Z, Y, X)) – Input 3-D volume.

Returns:

Radially averaged values. N = max(Z, Y, X) // 2 + 1.

Return type:

ndarray, shape (N,)

susan.utils._functions.fsc_get(v1, v2, msk=None)

Compute the Fourier Shell Correlation (FSC) between two half-maps.

\[FSC(r) = \frac{ \text{RadialAvg}_r\!\left( \mathcal{F}\{v_1 \cdot m\} \cdot \overline{\mathcal{F}\{v_2 \cdot m\}} \right) }{\sqrt{ \text{RadialAvg}_r\!\left(|\mathcal{F}\{v_1 \cdot m\}|^2\right) \cdot \text{RadialAvg}_r\!\left(|\mathcal{F}\{v_2 \cdot m\}|^2\right) }}\]

where m is the mask (1 everywhere if not provided).

Parameters:

• v1 (ndarray or str) – Input half-maps. Can be 3-D numpy arrays or paths to MRC files. Both must have the same shape.

• v2 (ndarray or str) – Input half-maps. Can be 3-D numpy arrays or paths to MRC files. Both must have the same shape.

• msk (ndarray or str or None, optional) – Real-space mask applied to both half-maps before the FFT. Can be a numpy array or a path to an MRC file. None (default) uses no mask.

Returns:

FSC curve. Shell 0 is set to 1.0; N = v1.shape[2] // 2 + 1.

Return type:

ndarray, shape (N,)

susan.utils._functions.fsc_analyse(fsc, apix=1.0, thres=0.143)

Find the resolution where the FSC drops below a threshold.

Parameters:

• fsc (array_like) – FSC curve as returned by fsc_get.

• apix (float or array_like, optional) – Pixel size in Angstroms. Default 1.0 (returns resolution in pixels).

• thres (float, optional) – FSC threshold. Default 0.143 (gold-standard half-map criterion).

Returns:

Named tuple with fields:

• fpix — resolution in Fourier pixels (int).

• res — resolution in Angstroms (float). 0.0 if the FSC never drops below thres.

Return type:

datatypes.fsc_info

susan.utils._functions.fsc_get_fpix(fsc, th_list=(0.5, 0.143), interp=True)

Resolution in Fourier pixels at each of several FSC thresholds.

Unlike fsc_analyse(), this always returns a list, one entry per threshold, and can interpolate the crossing to sub-shell precision.

Parameters:

• fsc (array_like) – FSC curve as returned by fsc_get() (n = box//2 + 1 shells).

• th_list (float or sequence of float, optional) – FSC threshold(s). A scalar is accepted and treated as a 1-element list. Default (0.5, 0.143).

• interp (bool, optional) – If True (default) linearly interpolate between the two shells bracketing the crossing. Matters when the crossings are only a few shells apart, which is typical of CryoET half-map FSCs: with integer crossings the anchors used by ssnr_from_fsc() can be off by more than 15%.

Returns:

One entry per threshold, always a list even for a single threshold.

Sentinels:

• nan — the FSC never drops below the threshold.

• 0.0 — the FSC is already below the threshold at shell 0.

Return type:

list of float

susan.utils._functions.ssnr_from_fsc(fsc, apix, th_list=(0.5, 0.143), n_eff=None, fallback=True)

Estimate the ad-hoc SSNR parameters (S, F) from an FSC curve.

SUSAN models the spectral SNR as (see susan.utils.datatypes.ssnr)

\[SSNR(s) = 10^{3S} \cdot e^{-100 F s}, \quad s \text{ in } 1/\text{\AA}\]

so \(\ln SSNR\) is linear in s and two points determine it. The two anchors are taken from the FSC itself, converting each threshold t to the map SSNR it corresponds to, \(SSNR = 2t/(1-t)\):

\[F = \frac{\ln(q_1/q_2)\,(box \cdot apix)}{100\,(r_2-r_1)}, \quad S = \frac{\ln q_1 + \ln(q_1/q_2)\, r_1/(r_2-r_1)}{3 \ln 10}\]

with \(r_i\) the crossing radii in Fourier pixels and \(q_i = 2t_i/(1-t_i)\). Note that S depends only on the ratio \(r_1/(r_2-r_1)\) and so is independent of the pixel size; F scales with box*apix.

The two-anchor solve is used rather than a least-squares fit because the applied weight, SSNR/(1+SSNR), saturates at 0 and 1: only the location and sharpness of the turnover matter, and those are what the anchors fix.

Parameters:

• fsc (array_like) – FSC curve as returned by fsc_get() (n = box//2 + 1 shells, so the box side is 2*(n-1)).

• apix (float) – Pixel size in Angstroms of the maps the FSC was computed from. Needed for F; S does not depend on it.

• th_list (sequence of two float, optional) – The two FSC anchors, in decreasing order. Default (0.5, 0.143).

• n_eff (float or None, optional) –

  Effective number of independent 2D measurements contributing to the map, roughly n_particles * n_tilts (times the symmetry order; use sum(prj_w) in place of n_tilts if the weights are not all 1). The FSC measures the SSNR of the map, while the model is defined as the SSNR of a single projection, so S is reduced by log10(n_eff)/3 (F is unchanged). Pass it when the result is destined for the substack whitening or the reconstruction Wiener filter, both of which apply the SSNR once per projection. None (default) returns the map SSNR unconverted.

  Note that with n_eff=None the returned S is always positive (every term of its formula is, as long as th_list[0] > 1/3), so the 10^(3S) > 1 gate in the C++ radial_frc_acc always engages. Applying the conversion can push S below 0 on a low-SSNR dataset, which that gate reads as “disabled” and silently reverts to pure whitening.

• fallback (bool, optional) –

  What to do when the requested anchors are not both reached. Default True, which walks down this ladder:

  1. Both th_list crossings usable — solve directly.

  2. Only th_list[0] crossed — re-anchor on ((1+th_list[0])/2, th_list[0]), i.e. (0.75, 0.5) for the defaults. The estimate is then extrapolated past the data, and a warning is issued.

  3. Neither crossed — the map is Nyquist-limited rather than SNR-limited. Anchor on the curve itself (first shell below 0.99 and the outermost shell) so the taper follows the FSC value at Nyquist. It is correspondingly mild, tending to no taper at all as the FSC flattens. A warning is issued.

  Set False to get nan instead of any fallback.

Returns:

Named tuple with fields S and F.

Both fields are nan if no usable pair of anchors exists: the crossings are out of order or coincident, the first anchor sits at shell 0, the curve has fewer than 2 shells, or fallback is False and the requested anchors were not both reached. Check with math.isnan(rslt.F) before use.

Return type:

datatypes.ssnr

Raises:

ValueError – If th_list does not hold exactly two decreasing values in (0, 1), or if apix is not positive.

susan.utils._functions.bandpass(v, lowpass, highpass=0, rolloff=1)

Apply a bandpass filter to a 3-D volume in Fourier space.

Both the low-pass and high-pass edges use a cosine rolloff, giving a smooth (Hann-like) transition rather than a hard cut.

Parameters:

• v (ndarray, shape (Z, Y, X)) – Input volume.

• lowpass (float) – Low-pass cutoff in Fourier pixels (0 = no low-pass). Shells above this radius are attenuated.

• highpass (float, optional) – High-pass cutoff in Fourier pixels (0 = no high-pass, default). Shells below this radius are attenuated.

• rolloff (int, optional) – Width of the cosine rolloff in Fourier pixels. Default 1.

Returns:

Filtered volume, same shape as v.

Return type:

ndarray, float32

susan.utils._functions.apply_FOM(v, fsc_array)

Apply a figure-of-merit (FOM) filter derived from an FSC curve.

Multiplies each Fourier shell by \(\sqrt{FSC}\):

\[v_{\text{FOM}} = \mathcal{F}^{-1}\!\left\{ \mathcal{F}\{v\} \cdot \sqrt{\text{fsc\_array}} \right\}\]

See Rosenthal & Henderson (2003).

Parameters:

• v (ndarray, shape (Z, Y, X)) – Input volume.

• fsc_array (array_like, shape (N,)) – FSC curve as returned by fsc_get. Values are clipped to [0, 1] before taking the square root.

Returns:

FOM-weighted volume, same shape as v.

Return type:

ndarray, float32

susan.utils._functions.spectral_weight_CFSC(v, apix=0.0, ssnr=(0.0, 0.0))

Compute the CFSC radial spectral weight of a volume.

This is the Python equivalent of the cfsc whitening the GPU code applies to a 3-D reference (RadialAverager::preset_FRC_vol), returned as a profile instead of being applied in place, so the weight measured on one map can be applied to another with apply_spectral_weight().

For each shell r = round(|k|) of the half-spectrum:

\[w[r] = \sqrt{\sum_{|k| \in r} \frac{|F(k)|^2}{\max(r,1)}} \cdot \mathrm{gain}[r] \cdot \sqrt{N/2}\]

with w[0] = sqrt(N/2). The division by r makes the 3-D shell sum scale like the 2-D ring sum used on the substack, keeping the reference and the data whitened consistently. The optional ad-hoc SSNR gain is (s+1)/max(s,1e-6) with s = 10^(3 S) exp(-100 F r / (N apix)), active only when 10^(3 S) > 1.

No masking is performed: apply the mask to v beforehand if the weight should be measured over the masked region only.

Parameters:

• v (ndarray, shape (N, N, N)) – Input volume. Must be cubic with an even box size.

• apix (float, optional) – Pixel size in Angstroms, used only by the SSNR gain. Default: 0 (gain disabled).

• ssnr (tuple of float, optional) – Ad-hoc SSNR parameters (F, S), matching -ssnr_param. Default: (0, 0) (gain disabled).

Returns:

Radial weight to be divided out, as consumed by apply_spectral_weight().

Return type:

ndarray, float32, shape (N//2+1,)

susan.utils._functions.apply_spectral_weight(v, wgt)

Divide a volume by a radial spectral weight in Fourier space.

The counterpart of spectral_weight_CFSC(). Note the direction: the weight is divided out, unlike apply_FOM() and bandpass(), which multiply. Shells where the weight is not positive are zeroed, as is the DC term and everything beyond the end of wgt, matching radial_frc_norm_vol. The result therefore has zero mean.

No masking is performed: re-apply the mask afterwards if needed.

Parameters:

• v (ndarray, shape (N, N, N)) – Volume to weight. Must be cubic with an even box size.

• wgt (array_like, shape (N//2+1,)) – Radial weight to divide out, typically from spectral_weight_CFSC().

Returns:

Weighted volume, same shape as v.

Return type:

ndarray, float32

susan.utils._functions.phase_randomize(v, fpix=0, seed=None)

Randomize the Fourier phases of a volume, preserving its amplitudes.

Builds a null map that shares the radial (and full 3-D) amplitude spectrum of v but carries no structural information above fpix. Correlating data against such a map measures the correlation obtainable by chance, which is the empirical noise floor for CC-based resolution cut-offs.

Phases are borrowed from the transform of a random real volume rather than drawn directly, so Hermitian symmetry is satisfied by construction and the result is exactly real. Because the amplitudes are untouched, spectral weighting and energy normalisation behave identically on the randomized map and on the original.

No masking is performed. Apply the same mask to the result that the original carries if the comparison is meant to be like for like.

Parameters:

• v (ndarray, shape (N, N, N)) – Input volume. Must be cubic with an even box size.

• fpix (float, optional) – Radius in Fourier pixels at which randomization starts. Shells with r >= fpix are randomized, shells below keep their original phases. Default: 0, which randomizes every frequency including DC.

• seed (int or None, optional) – Seed for the random phases. Default: None (non-reproducible).

Returns:

Phase-randomized volume, same shape as v.

Return type:

ndarray, float32

susan.utils._functions.fsc_sharpen(v, fsc, apix, bfactor, fom='rosenthal', lowpass=True, thres=0.143, rolloff=2)

FSC-weighted B-factor sharpening of a map from its half-map FSC.

Builds a single radial Fourier weight combining up to three per-shell terms and applies it to v:

1. B-factor amplification exp(-bfactor * s^2 / 4) where s = k / (N * apix) is the spatial frequency (1/Angstrom) of shell k. A negative bfactor sharpens (boosts high frequencies); a positive one blurs.

2. FOM weighting derived from the FSC, which tapers the amplification to zero where the half-maps stop correlating, so noise beyond the resolution limit is not amplified. 'rosenthal' uses the full-map figure of merit \(\sqrt{2\,FSC/(1+FSC)}\) (appropriate for the combined map; Rosenthal & Henderson, 2003); 'sqrt' uses \(\sqrt{FSC}\) (matches apply_FOM()); None disables it.

3. Cosine low-pass at the FSC resolution (fsc_analyse() with thres), a safety cut so shells past the resolution are not boosted.

This is the principled alternative to a blind (FSC-agnostic) B-factor: the amplification is gated by where there is real, reproducible signal.

Parameters:

• v (ndarray, shape (Z, Y, X)) – Map to sharpen (typically the combined reconstruction).

• fsc (array_like, shape (N//2+1,)) – Half-map FSC curve, as returned by fsc_get().

• apix (float) – Pixel size in Angstroms.

• bfactor (float) – B-factor in Angstrom^2. Negative sharpens, positive blurs.

• fom ({'rosenthal', 'sqrt', None}, optional) – FSC figure-of-merit weighting. Default 'rosenthal'.

• lowpass (bool, optional) – Apply a cosine low-pass at the FSC resolution. Default True.

• thres (float, optional) – FSC threshold used to locate the low-pass edge. Default 0.143.

• rolloff (int, optional) – Cosine rolloff width (Fourier pixels) of the low-pass. Default 2.

Returns:

Sharpened volume, same shape as v.

Return type:

ndarray, float32

See also

apply_FOM

FSC figure-of-merit weighting only (no B-factor).

fsc_sharpen_filter

wrap this as a map_filter_fsc callable for SubtomoAvg.

susan.utils._functions.fsc_sharpen_filter(apix, bfactor, **kwargs)

Build a map_filter_fsc callable that applies fsc_sharpen().

The returned function has signature filter(vol, fsc) -> vol, matching SubtomoAvg.map_filter_fsc, so it plugs directly into the post-processing hook — the project passes the per-reference FSC in automatically each iteration.

Parameters:

• apix (float) – Pixel size in Angstroms (e.g. sta.pix_size).

• bfactor (float) – B-factor in Angstrom^2 (negative sharpens).

• **kwargs – Forwarded to fsc_sharpen() (fom, lowpass, thres, rolloff).

Returns:

lambda vol, fsc: fsc_sharpen(vol, fsc, apix, bfactor, **kwargs).

Return type:

callable

Examples

>>> sta.map_filter_fsc = fsc_sharpen_filter(sta.pix_size, bfactor=-120)

susan.utils._functions.get_extension(filename)

Return the file extension including the leading dot.

Parameters:

filename (str)

Returns:

Extension, e.g. '.mrc'. Empty string if there is no extension.

Return type:

str

susan.utils._functions.is_extension(filename, extension)

Check whether filename has the given extension (case-sensitive).

Parameters:

• filename (str)

• extension (str) – With or without a leading dot (both forms are accepted).

Return type:

bool

susan.utils._functions.force_extension(filename, extension)

Return filename with its extension replaced by extension.

Parameters:

• filename (str)

• extension (str) – With or without a leading dot (both forms are accepted).

Returns:

Path with the new extension.

Return type:

str

susan.utils._functions.time_now()

Return the current local date and time.

Return type:

datetime.datetime

susan.utils._functions.create_sphere(r, N, center=None)

Create a soft spherical mask of radius r in a cube of side N.

The mask value at each voxel is clip(r - radius, 0, 1), giving a smooth 1-pixel-wide transition at the sphere boundary.

Parameters:

• r (float) – Sphere radius in pixels.

• N (int) – Side length of the output cube.

• center (array-like of 3 floats, optional) – Voxel coordinates (z, y, x) of the sphere centre. None (default) places it at the geometric centre (N//2, N//2, N//2), reproducing the original behaviour. Use this to place a mask at an arbitrary location (e.g. a tile centre in local processing).

Returns:

Soft spherical mask; 1 inside, 0 outside, linear transition at edge.

Return type:

ndarray, shape (N, N, N), float32

susan.utils._functions.bin_vol(vol, bin_level)

Low-pass filter and downsample a volume by a power of two.

Applies a low-pass filter at the new Nyquist frequency before downsampling to prevent aliasing.

Parameters:

• vol (ndarray, shape (N, N, N)) – Input volume.

• bin_level (int) – Downsampling factor as a power of two. bin_level=1 halves each dimension; bin_level=2 quarters it, etc.

Returns:

Downsampled volume of shape (N//s, N//s, N//s) where s = 2**bin_level.

Return type:

ndarray, float32

susan.utils._functions.bin_frame(in_frame, scale, out_frame=None)

Area-weighted downsample of a single 2-D frame by a float scale.

Output dimensions are ceil(H/scale) and ceil(W/scale). The window offset (N - N_b*scale)/2 - (scale-1)/2 keeps the sampling origin on SUSAN’s pixel-centre convention (input index i at coordinate i, tomogram centre at stk_center = N/2), so the same particle position projects to the same physical point across binning levels. Preserving the geometric box-edge centre instead would shift binned content by (scale-1)/2 input pixels and blur the reconstruction across tilts.

Edge bins extend past the input boundary; out-of-bounds contributions are skipped and each output pixel is normalised by the actual in-bounds weight, so no artificial padding is introduced.

Parameters:

• in_frame (ndarray, shape (H, W)) – Input frame; converted to contiguous float32 if needed.

• scale (float, > 1.0) – Downsampling factor (input pixels per output pixel).

• out_frame (ndarray, optional) – Pre-allocated output buffer of shape (ceil(H/scale), ceil(W/scale)), dtype float32, contiguous. Allocated internally if not given.

Returns:

The downsampled frame.

Return type:

ndarray, float32, shape (ceil(H/scale), ceil(W/scale))

susan.utils._functions.bin_frame_shape(H, W, scale)

Return the (H_b, W_b) shape that bin_frame() would produce.

susan.utils._functions.mask_diameter(mask_file, threshold=0.5)

Estimate the particle diameter in pixels from a soft mask MRC file.

The diameter is that of the sphere whose volume equals the volume of mask voxels above threshold. Returning pixels (not Angstroms).

Parameters:

• mask_file (str) – Path to the mask MRC file.

• threshold (float, optional) – Voxel values above this level are considered ‘inside’ the mask. Default 0.5 works for all standard soft masks.

Returns:

Equivalent-sphere diameter in pixels.

Return type:

float

susan.utils._functions.angular_step_from_fsc(fsc_fpix)

Angular step from an FSC resolution in Fourier pixels.

Returns the angle subtended by one Fourier pixel at the resolution shell fsc_fpix:

Δθ = atan2(1, fsc_fpix)   [degrees]

This is the smallest orientation change that moves the projected signal by one pixel at the resolution limit — i.e. the Nyquist angular step for the given resolution. No pixel size or particle diameter is needed.

Parameters:

fsc_fpix (int or float) – Resolution in Fourier pixels as returned by fsc_analyse.

Returns:

Suggested angular step in degrees.

Return type:

float

susan.utils._functions.is_odd(v)

Check whether an integer is odd.

Parameters:

v (int)

Returns:

True if v is odd, False otherwise.

Return type:

bool

susan.utils._functions.is_even(v)

Check whether an integer is even.

Parameters:

v (int)

Returns:

True if v is even, False otherwise.

Return type:

bool