diff --git a/README.md b/README.md index 6feb8df..a7c2e3d 100644 --- a/README.md +++ b/README.md @@ -56,8 +56,9 @@ Option | Default | Description Option | Default | Description --|--|-- `-m`, `--mode` | `angular` | Destriping approach: `angular` (recommended) or `linear` (deprecated). See [Destriping mode](#destriping-mode) below. -`--angle` | *(auto-estimated)* | Curtaining orientation, degrees from horizontal. Omit to estimate automatically from the tilt series' central frame[^estimation]; pass a value to override. A diagnostic plot is saved alongside the output when auto-estimated[^diagnosticplot]. -`--reference-frame` | *(mid-stack index)* | Stack index used for angle estimation and the destriping preview. See [Choosing a reference frame](#choosing-a-reference-frame) for more information. +`--angle` | *(auto-estimated)* | Curtaining orientation, degrees from horizontal. Omit to estimate automatically per frame[^estimation] (see [Per-frame curtain angle](#per-frame-curtain-angle) below); pass a value to override and skip estimation entirely. A diagnostic plot for the reference frame is saved alongside the output when auto-estimated[^diagnosticplot]. +`--reference-frame` | *(mid-stack index)* | Stack index used to seed the per-frame angle walk and as the destriping preview frame. See [Choosing a reference frame](#choosing-a-reference-frame) for more information. +`--anchor-tilts` | `5` | Number of frames nearest the reference frame used to seed the per-frame consensus walk. See [Per-frame curtain angle](#per-frame-curtain-angle) below. `--angular-width` | `8.0` | Angular width of the destriping notch, in degrees. Only used when `--mode angular`. Narrower keeps more real structure sharing a nearby angle to the curtains, at the cost of weaker curtain removal. `--notch-fraction` | `0.03` | Width of the destriping notch, as a fraction of image width. Narrower removes less real signal running parallel to the curtains, but leaves more curtaining behind. Only used when `--mode linear`, which is deprecated. `--protect-fraction` | `0.01` | Fraction of image width around the zero-frequency (DC) origin exempted from destriping. See [Destriping mode](#destriping-mode) below for why this exists and its trade-off. Only used when `--mode linear`, which is deprecated. @@ -66,7 +67,7 @@ Option | Default | Description Option | Default | Description --|--|-- `--output-dir` | *(required for directory input)* | Output directory for batch mode, mirroring the input directory's structure. -`--angle-outlier-threshold` | `5.0` | Warn if an individual series' own angle estimate differs from the batch consensus by more than this many degrees. +`--angle-outlier-threshold` | `5.0` | Within a stack, the maximum degrees a frame's own angle estimate may deviate from its nearest resolved neighbor before falling back to that neighbor's angle (see [Per-frame curtain angle](#per-frame-curtain-angle)). In batch mode, also used to warn (only) when a series' reference-frame estimate differs from the cross-stack diagnostic consensus by more than this many degrees. `--workers` | `0` (all CPUs) | Number of parallel processes to use in batch mode. ##### Other options @@ -110,19 +111,30 @@ pylisc stack raw_tilt_series/ --output-dir cleared_tilt_series/ Files already carrying a `_PyLisC_` suffix (i.e. previous PyLisC output) are skipped, so re-running against the same directory won't reprocess its own results. -#### Shared curtain angle -Unless `--angle` is given explicitly, single-file mode estimates the curtaining angle from one frame. In batch mode, PyLisC instead estimates an angle **per series** and combines them into a single shared angle, which is then applied to every series in the batch rather than letting each one drift independently. +#### Per-frame curtain angle +As curtaining angles can drift slightly across a tilt series, PyLisC only applies a single angle when one is provided explicitly with `--angle`. Otherwise, for each tilt series PyLisC: +1. Estimates each frame's curtaining angle. +2. Determines a 'seed' angle from the `--anchor-tilts` frames nearest the [reference frame](#choosing-a-reference-frame). +3. Walks outward from the seed frames, comparing each frame's estimate to the nearest already-resolved frame and either: + 1. Accepts the estimate: if it's within `--angle-outlier-threshold` degrees of the nearest already-resolved frame (itself becoming the new "nearest resolved" frame). + 2. Rejects the estimate: replacing with that nearest resolved angle and logged as a warning naming the frame. + - This allows for slight angle drift over a series while preventing one bad frame from throwing off its neighbours. +4. Destripes each frame using the resolved angle. -The combination is a confidence-weighted circular mean: each series' angle is weighted by its own confidence ratio (see [Curtain angle diagnostic plot](#curtain-angle-diagnostic-plot) below), so a series with a clear, sharp peak counts for more than one with a flat, uncertain profile. +The reference frame's diagnostic plot (see [Curtain angle diagnostic plot](#curtain-angle-diagnostic-plot) below) is saved alongside the output. -#### Outlier detection +#### Multi-series consensus check -If any individual series' own angle estimate differs from the batch consensus by more than `--angle-outlier-threshold` (default `5.0` degrees), a warning is printed naming that series: +In batch mode, PyLisC additionally combines each tilt series' reference-frame estimate into a consensus angle to use as a diagnostic value. This is a confidence-weighted circular mean: each series' angle is weighted by its own confidence ratio (see [Curtain angle diagnostic plot](#curtain-angle-diagnostic-plot) below), clipped to a per-run cap first so a single sharp FFT peak can't dominate the consensus (the cap is the confidence distribution's median plus 3x its scaled median absolute deviation, falling back to 5x the median when every value is identical). + +If any individual series' own reference-frame estimate differs from this consensus by more than `--angle-outlier-threshold` (default `5.0` degrees), a warning is printed naming that series: ``` WARNING: sample_07.mrc angle (58.3 deg) deviates 41.2 deg from consensus (17.1 deg) -- check its diagnostic plot ``` +This does not affect the actual destriping, however, with each tilt series destriped using its own per-frame walk described [above](#per-frame-curtain-angle). + ### `pylisc frames` ```sh @@ -132,11 +144,10 @@ pylisc frames [OPTIONS] --output-dir OUTPUT_DIR --filename-template TEMPLATE INP This command is aimed at destriping tilt images that exist as individual 2D MRC frames, i.e. not yet assembled/aligned into a stack. `INPUT_DIR` is not recursed into; every `*.mrc` directly inside it (excluding PyLisC's own `_PyLisC_` output) is treated as one tilt image. `--output-dir` is required, and mirrors the input's flat structure: each frame is written back out individually with a `_PyLisC_{mode}` suffix, same as `--output-dir` does for [batch mode](#batch-mode). #### Options -`pylisc frames` uses many of the same options as `pylisc stack`, see [above](#options) or run `pylisc frames -h` for further information. Frames mode does not take `--reference-frame` or `--preview-strengths`, but adds: +`pylisc frames` uses many of the same options as `pylisc stack`, see [above](#options) or run `pylisc frames -h` for further information, including `--anchor-tilts` (here, the number of tilt buckets rather than frames used to seed the consensus walk — see [Per-tilt curtain angle](#per-tilt-curtain-angle) below). Frames mode does not take `--reference-frame` or `--preview-strengths`, but adds: Option | Default | Description --|--|-- -`--anchor-tilts` | `5` | Number of tilt buckets nearest each series' median tilt used to seed its consensus walk. See [Per-tilt curtain angle](#per-tilt-curtain-angle) below. `--print-angles` | off | Print a diagnostic table of per-file and per-bucket angle estimation, and write it to CSV files in `--output-dir`. #### Filename template @@ -158,7 +169,7 @@ Field boundaries default to underscore only. `--filename-delimiters` sets which Curtaining orientation drifts slightly with tilt angle, so unless `--angle` is given explicitly, frames mode does **not** use one consensus angle for the whole directory. Instead: 1. Every frame's own angle is estimated. 2. Each frame's confidence ratio is clipped to a per-run cap before use, so a single sharp FFT peak can't dominate its bucket's consensus. The cap is the confidence distribution's median plus 3x its (scaled) median absolute deviation, falling back to 5x the median when every value is identical, (i.e. zero deviation). -3. Frames are grouped by tilt angle, rounded to the nearest whole degree (so e.g. two positions' `-30.00°` and `-29.98°` tilts fall in the same bucket), and each bucket's estimates are combined into a per-tilt consensus (the same confidence-weighted circular mean as [batch mode](#shared-curtain-angle)), using the clipped confidences. +3. Frames are grouped by tilt angle, rounded to the nearest whole degree (so e.g. two positions' `-30.00°` and `-29.98°` tilts fall in the same bucket), and each bucket's estimates are combined into a per-tilt consensus (the same confidence-weighted circular mean used in [stack mode's cross-stack consensus check](#multi-series-consensus-check)), using the clipped confidences. 4. Tilt buckets are clustered into separate acquisition series (by matching frame count per bucket, then splitting wherever the tilt spacing breaks step), so a directory holding more than one tilt series is handled independently per series rather than as one pool. 5. Within each series, the `--anchor-tilts` buckets around the series' median tilt are used as a trusted seed. PyLisC then walks outward from this seed, and compares each tilt bucket's own consensus to the nearest already-resolved bucket. A tilt bucket's consensus angle is accepted (becoming the new "nearest resolved" point) if it's within `--angle-outlier-threshold` degrees of it, otherwise it is replaced with that nearest resolved angle and logged as a warning naming both tilts. This approach prevents buckets with incorrect angles from throwing off subsequent tilts, while allowing for slight angle drift over a tilt series. 6. Pass `--print-angles` to see the full per-file and per-bucket breakdown (angle, confidence, resolved angle, seed/accepted/rejected status) as tables and CSV files written to `--output-dir`. @@ -193,7 +204,7 @@ Curtaining removal works by finding curtaining's signature in Fourier space and - **`linear` (deprecated).** Dims frequencies by their *distance* from the curtain line rather than their direction. Below a radius set by `--notch-fraction`, distance alone can no longer distinguish direction at all, so without `--protect-fraction` exempting a small disc around the origin, large-scale contrast gets suppressed at every angle near that radius, not just along the curtains. Protecting that disc, in turn, risks letting broad, low-frequency curtaining pass through unfiltered if the curtaining's own frequency sits close to the protected radius. `angular` avoids this trade-off entirely. ### Choosing a reference frame -The reference frame should be the tilt with the least foreshortening and the best signal-to-noise, since that gives the most reliable curtain angle estimate and the clearest destriping preview. In practice this is the 0° tilt (or the pretilt used during lamella imaging). +The reference frame should be the tilt with the least foreshortening and the best signal-to-noise, since that gives the most reliable curtain angle estimate, the clearest destriping preview, and the most trustworthy seed for the [per-frame consensus walk](#per-frame-curtain-angle). In practice this is the pre-tilt used during lamella imaging. - **Dose-symmetric schemes** (0° acquired first, then alternating ±): use `--reference-frame 0`. - **Continuous sweeps** (most-negative tilt acquired first): 0° sits in the middle of the stack, so use roughly `--reference-frame ` for an n-tilt series. diff --git a/pyproject.toml b/pyproject.toml index 2f248ed..6250daa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pylisc" -version = "2.4.0" +version = "2.5.0" description = "Python implementation of LisC algorithm" readme = "README.md" authors = [ diff --git a/src/pylisc/cli.py b/src/pylisc/cli.py index bc6c621..bdb5001 100644 --- a/src/pylisc/cli.py +++ b/src/pylisc/cli.py @@ -91,7 +91,7 @@ def version_callback(value: bool | None) -> None: ] AnchorTiltsOpt = Annotated[ int, - typer.Option('--anchor-tilts', help='Number of tilt buckets nearest each series\' median tilt used to seed its consensus walk (frames mode only).', rich_help_panel='Batch options', min=1), + typer.Option('--anchor-tilts', help='Number of tilt buckets (frames mode) or frames (stack mode) nearest the median/reference used to seed the angle consensus walk.', rich_help_panel='Batch options', min=1), ] # Define callback for pylisc (to allow version option) @@ -139,6 +139,7 @@ def stack( notch_frac: NotchFracOpt = 0.03, dc_protect_frac: DcProtectFracOpt = 0.01, angle_outlier_threshold: AngleOutlierThresholdOpt = 5.0, + anchor_tilts: AnchorTiltsOpt = 5, version: VersionOpt = None, workers: WorkersOpt = 0, ): @@ -181,6 +182,7 @@ def stack( notch_frac=notch_frac, dc_protect_frac=dc_protect_frac, angle_outlier_threshold=angle_outlier_threshold, + anchor_tilts=anchor_tilts, preview_strengths=preview_strengths, force=force, dry_run=dry_run, diff --git a/src/pylisc/estimate_angle.py b/src/pylisc/estimate_angle.py index 42fc843..0ee1c26 100644 --- a/src/pylisc/estimate_angle.py +++ b/src/pylisc/estimate_angle.py @@ -89,6 +89,54 @@ def plot_angular_energy( plt.close(fig) return output_path +def clip_confidence_outliers(confidences: dict) -> dict: + ''' + Cap confidences above a median+MAD threshold so a single spuriously sharp FFT peak can't dominate a consensus + ''' + conf_values = np.array(list(confidences.values())) + if len(conf_values) == 0: + return dict(confidences) + median_conf = np.median(conf_values) + mad = np.median(np.abs(conf_values - median_conf)) + confidence_cap = median_conf + 3 * 1.4826 * mad if mad > 0 else median_conf * 5 + if confidence_cap <= 0: + return dict(confidences) + return {k: min(v, confidence_cap) for k, v in confidences.items()} + +def resolve_walk(keys_sorted: list, values: dict, angle_outlier_threshold: float, anchor_window: int) -> tuple[dict, dict]: + ''' + Seed a window of keys around the median key, and walk outward, checking each key's own value against the nearest already-resolved (seed/accepted) value + ''' + if len(keys_sorted) == 1: + k = keys_sorted[0] + return {k: values[k]}, {k: 'seed'} + median_key = np.median(keys_sorted) + center_idx = int(np.argmin([abs(k - median_key) for k in keys_sorted])) + window = max(1, anchor_window) + half = window // 2 + start = max(0, center_idx - half) + end = min(len(keys_sorted), start + window) + start = max(0, end - window) + seed_keys = keys_sorted[start:end] + resolved = {k: values[k] for k in seed_keys} + status = {k: 'seed' for k in seed_keys} + for direction, idx, edge in ((-1, start - 1, start), (1, end, end - 1)): + nearest = resolved[keys_sorted[edge]] + i = idx + while 0 <= i < len(keys_sorted): + k = keys_sorted[i] + own_value = values[k] + deviation = min(abs(own_value - nearest), 180 - abs(own_value - nearest)) + if deviation <= angle_outlier_threshold: + resolved[k] = own_value + status[k] = 'accepted' + nearest = own_value + else: + resolved[k] = nearest + status[k] = 'rejected' + i += direction + return resolved, status + def combine_angles(angles_deg: list, confidences: list) -> tuple: ''' Confidence-weighted circular mean of curtain angles diff --git a/src/pylisc/frames.py b/src/pylisc/frames.py index 692f609..51208a4 100644 --- a/src/pylisc/frames.py +++ b/src/pylisc/frames.py @@ -9,7 +9,7 @@ from rich.table import Table # Import internal PyLisC modules -from pylisc.estimate_angle import combine_angles, estimate_curtain_angle +from pylisc.estimate_angle import clip_confidence_outliers, combine_angles, estimate_curtain_angle, resolve_walk from pylisc.io import find_input_files, readMrcFile, writeMrcFile from pylisc.lisc import lisc_clear_frame from pylisc.log import logger, per_file_log @@ -178,42 +178,14 @@ def _resolve_series(series_buckets, bucket_consensus, angle_outlier_threshold, a Seed buckets around the median tilt, and walk out, checking each against the nearest already-resolved (seed/accepted) bucket ''' buckets_sorted = sorted(series_buckets) - if len(buckets_sorted) == 1: - b = buckets_sorted[0] - return {b: bucket_consensus[b]}, {b: 'seed'} - - median_tilt = np.median(buckets_sorted) - center_idx = int(np.argmin([abs(b - median_tilt) for b in buckets_sorted])) - window = max(1, anchor_tilts) - half = window // 2 - start = max(0, center_idx - half) - end = min(len(buckets_sorted), start + window) - start = max(0, end - window) - - seed_buckets = buckets_sorted[start:end] - resolved = {b: bucket_consensus[b] for b in seed_buckets} - status = {b: 'seed' for b in seed_buckets} - - for direction, idx, edge in ((-1, start - 1, start), (1, end, end - 1)): - nearest = resolved[buckets_sorted[edge]] - i = idx - while 0 <= i < len(buckets_sorted): - b = buckets_sorted[i] - own_angle = bucket_consensus[b] - deviation = min(abs(own_angle - nearest), 180 - abs(own_angle - nearest)) - if deviation <= angle_outlier_threshold: - resolved[b] = own_angle - status[b] = 'accepted' - nearest = own_angle - else: - resolved[b] = nearest - status[b] = 'rejected' - logger.warning( - "tilt {}° consensus angle ({}°) deviates {}° from nearest resolved angle ({}°) - using that angle instead", - b, f'{own_angle:.1f}', f'{deviation:.1f}', f'{nearest:.1f}', - ) - i += direction - + resolved, status = resolve_walk(buckets_sorted, bucket_consensus, angle_outlier_threshold, anchor_tilts) + for b in buckets_sorted: + if status[b] == 'rejected': + deviation = min(abs(bucket_consensus[b] - resolved[b]), 180 - abs(bucket_consensus[b] - resolved[b])) + logger.warning( + "tilt {}° consensus angle ({}°) deviates {}° from nearest resolved angle ({}°) - using that angle instead", + b, f'{bucket_consensus[b]:.1f}', f'{deviation:.1f}', f'{resolved[b]:.1f}', + ) return resolved, status def _estimate_per_tilt_angles(paths, tilt_of, angle_outlier_threshold, anchor_tilts=5, print_angles=False, output_dir=None): @@ -230,18 +202,10 @@ def _estimate_per_tilt_angles(paths, tilt_of, angle_outlier_threshold, anchor_ti # A single spuriously sharp FFT peak can otherwise dominate its bucket's consensus raw_confidences = dict(confidences) - conf_values = np.array(list(confidences.values())) - if len(conf_values): - median_conf = np.median(conf_values) - mad = np.median(np.abs(conf_values - median_conf)) - confidence_cap = median_conf + 3 * 1.4826 * mad if mad > 0 else median_conf * 5 - else: - confidence_cap = 0.0 - if confidence_cap > 0: - n_clipped = sum(1 for c in confidences.values() if c > confidence_cap) - if n_clipped: - logger.debug('clipping {} frame(s) with confidence above {}', n_clipped, f'{confidence_cap:.2f}') - confidences = {p: min(c, confidence_cap) for p, c in confidences.items()} + confidences = clip_confidence_outliers(confidences) + n_clipped = sum(1 for p in confidences if confidences[p] < raw_confidences[p]) + if n_clipped: + logger.debug('clipping {} frame(s) with confidence above cap', n_clipped) tilt_buckets = {} for path in paths: diff --git a/src/pylisc/stack.py b/src/pylisc/stack.py index 2d81b1e..7daf83e 100644 --- a/src/pylisc/stack.py +++ b/src/pylisc/stack.py @@ -8,7 +8,7 @@ from pathlib import Path # Import internal PyLisC modules -from pylisc.estimate_angle import combine_angles, estimate_curtain_angle, plot_angular_energy +from pylisc.estimate_angle import clip_confidence_outliers, combine_angles, estimate_curtain_angle, plot_angular_energy, resolve_walk from pylisc.io import find_input_files, readMrcFile, writeMrcFile from pylisc.lisc import lisc_clear_frame from pylisc.log import logger, per_file_log @@ -40,6 +40,8 @@ def _process_series( angular_width, notch_frac, dc_protect_frac, + angle_outlier_threshold, + anchor_tilts, force, dry_run, preview_strengths=None, @@ -66,14 +68,14 @@ def _process_series( # Estimate curtaining angle if not provided if curtain_angle is None: - resolved_angle, angular_energy = estimate_curtain_angle(data[resolved_reference_frame]) - if not dry_run: - plot_angular_energy(angular_energy, resolved_angle, output_dir=out_path.parent) - median_energy = np.median(angular_energy) - confidence = angular_energy.max() / median_energy if median_energy > 0 else 0.0 - logger.info('({}) estimated curtaining angle: {}° (confidence: {})', path.name, resolved_angle, confidence) + per_frame_angle = _estimate_per_frame_angles( + path, data, resolved_reference_frame, angle_outlier_threshold, anchor_tilts, out_path.parent, dry_run, + ) + resolved_angle = per_frame_angle[resolved_reference_frame] + logger.info('({}) reference frame ({}) resolved curtaining angle: {}°', path.name, resolved_reference_frame, resolved_angle) else: resolved_angle = curtain_angle + per_frame_angle = None if dry_run: logger.info('({}) [dry-run] would write to {} (angle: {}°, pixel size: {})', path.name, out_path, resolved_angle, resolved_pixel_size) @@ -97,14 +99,14 @@ def _process_series( logger.info('Strength preview saved to {}', preview_path) return None - # Apply LisC to each frame + # Apply LisC to each frame, using each frame's own resolved angle when estimated cleared_stack = np.empty_like(data, dtype=np.float32) for i, frame in enumerate(data): cleared_stack[i] = lisc_clear_frame( frame, decurtaining_mode=mode, pixel_size_nm=resolved_pixel_size, - curtain_angle=resolved_angle, + curtain_angle=per_frame_angle[i] if per_frame_angle is not None else resolved_angle, apply_filter=apply_filter, filter_threshold_nm=filter_threshold, angular_width_deg=angular_width, @@ -117,6 +119,36 @@ def _process_series( return resolved_angle +def _estimate_per_frame_angles(path, data, resolved_reference_frame, angle_outlier_threshold, anchor_tilts, output_dir, dry_run): + ''' + Estimate a curtain angle for every frame in a stack, then walk out from the reference frame accepting each frame's own + estimate only if it's within angle_outlier_threshold of its nearest already-resolved neighbor (frame index is already + tilt-ordered in an MRC stack, so neighboring frames should carry similar angles) + ''' + angles = {} + ref_energy = None + for i, frame in enumerate(data): + angle, energy = estimate_curtain_angle(frame) + angles[i] = angle + if i == resolved_reference_frame: + ref_energy = energy + logger.debug('({}) frame {} est. angle: {}', path.name, i, angle) + + if not dry_run and ref_energy is not None: + plot_angular_energy(ref_energy, angles[resolved_reference_frame], output_dir=output_dir) + + frame_indices = list(range(len(data))) + resolved, status = resolve_walk(frame_indices, angles, angle_outlier_threshold, anchor_tilts) + for i in frame_indices: + if status[i] == 'rejected': + deviation = min(abs(angles[i] - resolved[i]), 180 - abs(angles[i] - resolved[i])) + logger.warning( + "({}) frame {} angle ({}°) deviates {}° from nearest resolved angle ({}°) - using that angle instead", + path.name, i, f'{angles[i]:.1f}', f'{deviation:.1f}', f'{resolved[i]:.1f}', + ) + return resolved + + def run_stack( input_path, output_mrc, @@ -131,6 +163,7 @@ def run_stack( notch_frac, dc_protect_frac, angle_outlier_threshold, + anchor_tilts, force, dry_run, workers, @@ -150,6 +183,7 @@ def run_stack( notch_frac=notch_frac, dc_protect_frac=dc_protect_frac, angle_outlier_threshold=angle_outlier_threshold, + anchor_tilts=anchor_tilts, force=force, dry_run=dry_run, workers=workers, @@ -168,6 +202,8 @@ def run_stack( angular_width=angular_width, notch_frac=notch_frac, dc_protect_frac=dc_protect_frac, + angle_outlier_threshold=angle_outlier_threshold, + anchor_tilts=anchor_tilts, preview_strengths=preview_strengths, force=force, dry_run=dry_run, @@ -187,6 +223,7 @@ def _run_stack_batch( notch_frac, dc_protect_frac, angle_outlier_threshold, + anchor_tilts, force, dry_run, workers, @@ -197,7 +234,9 @@ def _run_stack_batch( logger.info('{} files found in {}', len(series_paths), input_dir) if curtain_angle is None: - angles, confidences = [], [] + # Each stack is its own independent acquisition, so its own per-frame estimate is used for processing (below). + # This pass is diagnostic only: a reference-frame consensus across files, to warn on stacks that stand out. + angles, confidences = {}, {} for path in series_paths: data, _ = readMrcFile(path) if data is None: @@ -206,20 +245,21 @@ def _run_stack_batch( frame = data[frame_index] angle, energy = estimate_curtain_angle(frame) median_energy = np.median(energy) - confidence = energy.max() / median_energy if median_energy > 0 else 0.0 - angles.append(angle); confidences.append(confidence) - logger.debug('({}) est. angle: {} (conf.: {})', path.name, angle, confidence) + confidences[path] = energy.max() / median_energy if median_energy > 0 else 0.0 + angles[path] = angle + logger.debug('({}) est. angle: {} (conf.: {})', path.name, angle, confidences[path]) - consensus_angle, agreement = combine_angles(angles, confidences) + # A single spuriously sharp FFT peak can otherwise dominate the batch consensus + clipped_confidences = clip_confidence_outliers(confidences) + + consensus_angle, agreement = combine_angles(list(angles.values()), list(clipped_confidences.values())) logger.info('batch consensus angle: {}° (agreement: {})', f'{consensus_angle:.1f}', f'{agreement:.3f}') - for path, angle in zip(series_paths, angles): + for path, angle in angles.items(): deviation = min(abs(angle - consensus_angle), 180 - abs(angle - consensus_angle)) if deviation > angle_outlier_threshold: logger.warning('({}) est. angle ({}°) deviates {}° from consensus ({}°) - check diagnostic plot', path.name, f'{angle:.1f}', f'{deviation:.1f}', f'{consensus_angle:.1f}') - curtain_angle = consensus_angle - jobs = [] for path in series_paths: relative = path.relative_to(input_dir) @@ -246,6 +286,8 @@ def _run_stack_batch( angular_width=angular_width, notch_frac=notch_frac, dc_protect_frac=dc_protect_frac, + angle_outlier_threshold=angle_outlier_threshold, + anchor_tilts=anchor_tilts, force=force, dry_run=dry_run, ): path diff --git a/tests/integration/test_cli.py b/tests/integration/test_cli.py index 24b1884..ec41430 100644 --- a/tests/integration/test_cli.py +++ b/tests/integration/test_cli.py @@ -53,6 +53,14 @@ def test_reference_frame_at_series_length_is_rejected_cleanly(self, tmp_path, sy assert result.exit_code == 0 assert 'IndexError' not in result.output + def test_anchor_tilts_accepted_for_per_frame_walk(self, tmp_path, synthetic_tilt_series, write_synthetic_mrc): + input_path = tmp_path / 'series.mrc' + write_synthetic_mrc(input_path, synthetic_tilt_series(n_tilts=3, angle_deg=20)) + + result = runner.invoke(pylisc, ['stack', str(input_path), '--mode', 'angular', '--anchor-tilts', '2']) + assert result.exit_code == 0 + assert (tmp_path / 'series_PyLisC_angular.mrc').exists() + class TestCliBatch: def test_batch_run_with_outlier_warning(self, tmp_path, synthetic_tilt_series, write_synthetic_mrc): input_dir = tmp_path / 'raw' diff --git a/tests/unit/test_estimate_angle.py b/tests/unit/test_estimate_angle.py new file mode 100644 index 0000000..c6ce963 --- /dev/null +++ b/tests/unit/test_estimate_angle.py @@ -0,0 +1,46 @@ +''' +PyLisC: unit tests for shared consensus-robustness helpers (confidence clipping, anchor walk) +''' + +# Import external libraries +import pytest + +# Import internal functions +from pylisc.estimate_angle import clip_confidence_outliers, resolve_walk + +class TestClipConfidenceOutliers: + def test_no_outlier_leaves_confidences_unchanged(self): + confidences = {'a': 1.0, 'b': 1.1, 'c': 0.9} + assert clip_confidence_outliers(confidences) == confidences + + def test_spike_is_capped_others_untouched(self): + confidences = {'a': 1.0, 'b': 1.1, 'c': 0.9, 'spike': 500.0} + clipped = clip_confidence_outliers(confidences) + assert clipped['spike'] < confidences['spike'] + assert clipped['a'] == confidences['a'] + assert clipped['b'] == confidences['b'] + assert clipped['c'] == confidences['c'] + + def test_empty_dict(self): + assert clip_confidence_outliers({}) == {} + + +class TestResolveWalk: + def test_single_key_seeds_itself(self): + resolved, status = resolve_walk([0], {0: 12.0}, angle_outlier_threshold=5.0, anchor_window=5) + assert resolved == {0: 12.0} + assert status == {0: 'seed'} + + def test_smooth_series_all_accepted(self): + keys = list(range(7)) + values = {k: 10.0 + k for k in keys} + resolved, status = resolve_walk(keys, values, angle_outlier_threshold=5.0, anchor_window=3) + assert resolved == values + assert all(s in ('seed', 'accepted') for s in status.values()) + + def test_outlier_falls_back_to_nearest_resolved(self): + keys = list(range(7)) + values = {0: 50.0, 1: 50.0, 2: 50.0, 3: 50.0, 4: 50.0, 5: 50.0, 6: -10.0} + resolved, status = resolve_walk(keys, values, angle_outlier_threshold=5.0, anchor_window=3) + assert status[6] == 'rejected' + assert resolved[6] == pytest.approx(50.0) diff --git a/tests/unit/test_stack.py b/tests/unit/test_stack.py index 30b599a..050dff1 100644 --- a/tests/unit/test_stack.py +++ b/tests/unit/test_stack.py @@ -3,9 +3,10 @@ ''' # Import external libraries -import pytest +import numpy as np, pytest # Import internal functions +from pylisc.log import logger from pylisc.stack import _process_series, run_stack class TestProcessSeries: @@ -18,7 +19,8 @@ def test_explicit_pixel_size_overrides_header(self, tmp_path, synthetic_tilt_ser input_path, out_path, mode='angular', apply_filter=False, filter_threshold=5000.0, pixel_size=0.5, curtain_angle=20.0, reference_frame=0, angular_width=8.0, notch_frac=0.02, - dc_protect_frac=0.01, force=False, dry_run=False, + dc_protect_frac=0.01, angle_outlier_threshold=5.0, anchor_tilts=5, + force=False, dry_run=False, ) assert angle == pytest.approx(20.0) assert out_path.exists() @@ -33,7 +35,8 @@ def test_invalid_pixel_size_raises(self, tmp_path, synthetic_tilt_series, write_ input_path, out_path, mode='angular', apply_filter=False, filter_threshold=5000.0, pixel_size=0.0, curtain_angle=20.0, reference_frame=0, angular_width=8.0, notch_frac=0.02, - dc_protect_frac=0.01, force=False, dry_run=False, + dc_protect_frac=0.01, angle_outlier_threshold=5.0, anchor_tilts=5, + force=False, dry_run=False, ) def test_explicit_curtain_angle_skips_estimation(self, tmp_path, synthetic_tilt_series, write_synthetic_mrc): @@ -45,7 +48,8 @@ def test_explicit_curtain_angle_skips_estimation(self, tmp_path, synthetic_tilt_ input_path, out_path, mode='angular', apply_filter=False, filter_threshold=5000.0, pixel_size=0.34, curtain_angle=33.0, reference_frame=0, angular_width=8.0, notch_frac=0.02, - dc_protect_frac=0.01, force=False, dry_run=False, + dc_protect_frac=0.01, angle_outlier_threshold=5.0, anchor_tilts=5, + force=False, dry_run=False, ) assert angle == 33.0 @@ -62,7 +66,7 @@ def test_no_series_found_raises(self, tmp_path): mode='angular', apply_filter=False, filter_threshold=5000.0, pixel_size=None, curtain_angle=None, reference_frame=None, angular_width=8.0, notch_frac=0.02, dc_protect_frac=0.01, - angle_outlier_threshold=5.0, force=False, dry_run=False, workers=0, + angle_outlier_threshold=5.0, anchor_tilts=5, force=False, dry_run=False, workers=0, ) def test_batch_one_bad_file_does_not_abort_others(self, tmp_path, synthetic_tilt_series, write_synthetic_mrc): @@ -78,6 +82,85 @@ def test_batch_one_bad_file_does_not_abort_others(self, tmp_path, synthetic_tilt mode='angular', apply_filter=False, filter_threshold=5000.0, pixel_size=None, curtain_angle=None, reference_frame=None, angular_width=8.0, notch_frac=0.02, dc_protect_frac=0.01, - angle_outlier_threshold=5.0, force=False, dry_run=False, workers=1, + angle_outlier_threshold=5.0, anchor_tilts=5, force=False, dry_run=False, workers=1, ) assert (output_dir / 'good_PyLisC_angular.mrc').exists() + + def test_confidence_spike_does_not_skew_batch_consensus(self, tmp_path, synthetic_tilt_series, write_synthetic_mrc): + import re + input_dir = tmp_path / 'raw' + input_dir.mkdir() + output_dir = tmp_path / 'cleared' + # three ordinary low-tilt series striped at 20deg, with slightly different noise so their confidences aren't identical + write_synthetic_mrc(input_dir / 'low_a.mrc', synthetic_tilt_series(n_tilts=2, angle_deg=20, noise_std=5.0)) + write_synthetic_mrc(input_dir / 'low_b.mrc', synthetic_tilt_series(n_tilts=2, angle_deg=20, noise_std=6.0)) + write_synthetic_mrc(input_dir / 'low_c.mrc', synthetic_tilt_series(n_tilts=2, angle_deg=20, noise_std=4.5)) + # a single series with a much sharper (higher-confidence) peak at a wildly different angle + write_synthetic_mrc(input_dir / 'spike.mrc', synthetic_tilt_series(n_tilts=2, angle_deg=-60, amplitude=600.0, noise_std=1.0)) + + messages = [] + sink_id = logger.add(messages.append, level='INFO') + try: + run_stack( + input_path=input_dir, output_mrc=None, output_dir=output_dir, + mode='angular', apply_filter=False, filter_threshold=5000.0, + pixel_size=None, curtain_angle=None, reference_frame=None, + angular_width=8.0, notch_frac=0.02, dc_protect_frac=0.01, + angle_outlier_threshold=5.0, anchor_tilts=5, force=False, dry_run=False, workers=1, + ) + finally: + logger.remove(sink_id) + + match = next(re.search(r'batch consensus angle: (-?\d+\.\d+)°', str(m)) for m in messages if 'batch consensus angle' in str(m)) + consensus = float(match.group(1)) + # the low-tilt cluster (3 files) should still pull the consensus closer to itself than to the single spiky outlier + assert abs(consensus - 20.0) < abs(consensus - (-60.0)) + + +class TestEstimatePerFrameAngles: + def test_smooth_drift_uses_each_frame_own_estimate(self, tmp_path, synthetic_frame): + from pylisc.stack import _estimate_per_frame_angles + drift = [10, 12, 14, 16, 18] + stack = np.stack([synthetic_frame(angle_deg=a, seed=i) for i, a in enumerate(drift)]) + resolved = _estimate_per_frame_angles( + tmp_path / 'series.mrc', stack, resolved_reference_frame=2, + angle_outlier_threshold=5.0, anchor_tilts=3, output_dir=tmp_path, dry_run=True, + ) + for i, expected in enumerate(drift): + assert resolved[i] == pytest.approx(expected, abs=2.0) + + def test_outlier_frame_falls_back_to_nearest_resolved(self, tmp_path, synthetic_frame): + from pylisc.stack import _estimate_per_frame_angles + angles = [50, 50, 50, 50, 50, -40] + stack = np.stack([synthetic_frame(angle_deg=a, seed=i) for i, a in enumerate(angles)]) + + messages = [] + sink_id = logger.add(messages.append, level='WARNING') + try: + resolved = _estimate_per_frame_angles( + tmp_path / 'series.mrc', stack, resolved_reference_frame=2, + angle_outlier_threshold=5.0, anchor_tilts=3, output_dir=tmp_path, dry_run=True, + ) + finally: + logger.remove(sink_id) + + assert resolved[5] == pytest.approx(50.0, abs=2.0) + assert any('deviates' in str(m) and 'nearest resolved angle' in str(m) for m in messages) + + +class TestProcessSeriesPerFrame: + def test_estimation_applies_per_frame_not_one_shared_angle(self, tmp_path, synthetic_frame, write_synthetic_mrc): + input_path = tmp_path / 'series.mrc' + stack = np.stack([synthetic_frame(angle_deg=a, seed=i) for i, a in enumerate([10, 12, 14, 16, 18])]) + write_synthetic_mrc(input_path, stack) + out_path = tmp_path / 'out.mrc' + + angle = _process_series( + input_path, out_path, mode='angular', apply_filter=False, + filter_threshold=5000.0, pixel_size=0.34, curtain_angle=None, + reference_frame=2, angular_width=8.0, notch_frac=0.02, + dc_protect_frac=0.01, angle_outlier_threshold=5.0, anchor_tilts=3, + force=False, dry_run=False, + ) + assert angle == pytest.approx(14.0, abs=2.0) + assert out_path.exists() diff --git a/uv.lock b/uv.lock index 35979d1..01a0c80 100644 --- a/uv.lock +++ b/uv.lock @@ -390,7 +390,7 @@ wheels = [ [[package]] name = "pylisc" -version = "2.4.0" +version = "2.5.0" source = { editable = "." } dependencies = [ { name = "loguru" },