diff --git a/README.md b/README.md index 0f7e158..2234cd8 100644 --- a/README.md +++ b/README.md @@ -144,9 +144,10 @@ 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. 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 group). -3. Each group's estimates are combined into a per-tilt consensus (same confidence-weighted circular mean as [batch mode](#shared-curtain-angle)), which is the angle applied to every frame in that group. -4. All per-tilt consensus angles are then combined into an overall consensus, weighted both by confidence and by `cos(tilt)` (sample thickness grows ~1/cos(tilt) so high tilt angles are less reliable), so they count for less than well-sampled low-tilt groups rather than skewing the overall consensus by an equal vote. -5. Any per-tilt consensus that still deviates from the overall consensus by more than `--angle-outlier-threshold` is treated as unreliable, and will be destriped at the consensus angle of its nearest reliable tilt (by tilt-angle distance) instead, logging a warning naming both tilts. If every tilt ends up flagged, PyLisC falls back to the overall consensus for all of them. +3. Each frame's confidence ratio is clipped to a per-run cap before use, so a single sharp FFT peak can't dominate its group's consensus or the overall weighting below. The cap is the confidence distribution's median plus 3x its (scaled) median absolute deviation, which stays robust even with few frames (falling back to 5x the median when every value is identical, i.e. zero deviation). +4. Each group's estimates are combined into a per-tilt consensus (same confidence-weighted circular mean as [batch mode](#shared-curtain-angle)), using the clipped confidences, which is the angle applied to every frame in that group. +5. All per-tilt consensus angles are then combined into an overall consensus, weighted both by the group's typical (median) confidence and by `cos(tilt)` (sample thickness grows ~1/cos(tilt) so high tilt angles are less reliable), so they count for less than well-sampled low-tilt groups rather than skewing the overall consensus by an equal vote. +6. Any per-tilt consensus that still deviates from the overall consensus by more than `--angle-outlier-threshold` is treated as unreliable, and will be destriped at the consensus angle of its nearest reliable tilt (by tilt-angle distance) instead, logging a warning naming both tilts. If every tilt ends up flagged, PyLisC falls back to the overall consensus for all of them. #### Pixel size Individual frame MRCs frequently lack a reliable pixel size in their header, so frames mode does not fall back to it. Pixel size is only needed for the optional high-pass filter — if `--apply-filter` is set, `--pixel-size` must be given explicitly, or PyLisC exits with an error. diff --git a/pyproject.toml b/pyproject.toml index 64a1fc5..ca8fb5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pylisc" -version = "2.2.0" +version = "2.3.0" description = "Python implementation of LisC algorithm" readme = "README.md" authors = [ diff --git a/src/pylisc/frames.py b/src/pylisc/frames.py index 04bb0fe..31da8d5 100644 --- a/src/pylisc/frames.py +++ b/src/pylisc/frames.py @@ -141,6 +141,20 @@ def _estimate_per_tilt_angles(paths, tilt_of, angle_outlier_threshold): angles[path] = angle logger.debug('({}) tilt {}° est. angle: {} (conf.: {})', path.name, tilt_of[path], angle, confidences[path]) + # A single spuriously sharp FFT peak can otherwise dominate its bucket's consensus and the overall weighted average + 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()} + tilt_buckets = {} for path in paths: bucket = round(tilt_of[path]) @@ -156,7 +170,7 @@ def _estimate_per_tilt_angles(paths, tilt_of, angle_outlier_threshold): ) bucket_consensus[bucket] = consensus # High-tilt frames carry less signal (sample thickness grows ~1/cos(tilt)) so reduce weighting for overall consensus - bucket_weight[bucket] = sum(bucket_confidences) * np.cos(np.deg2rad(bucket)) + bucket_weight[bucket] = np.median(bucket_confidences) * np.cos(np.deg2rad(bucket)) logger.info('tilt {}°: consensus angle {}° (agreement: {}, n={})', bucket, f'{consensus:.1f}', f'{agreement:.3f}', len(bucket_paths)) overall_consensus, overall_agreement = combine_angles( diff --git a/tests/fixtures.py b/tests/fixtures.py index ec0045b..a62f155 100644 --- a/tests/fixtures.py +++ b/tests/fixtures.py @@ -36,8 +36,8 @@ def write_synthetic_mrc(path, stack, pixel_size_nm=3.4): mrc.voxel_size = pixel_size_nm * 10 # nm -> Angstrom -def write_synthetic_frame(path, angle_deg=0.0, pixel_size_nm=3.4, seed=0): +def write_synthetic_frame(path, angle_deg=0.0, pixel_size_nm=3.4, seed=0, **frame_kwargs): ''' Write a single synthetic 2D frame (as used in frames mode) ''' - write_synthetic_mrc(path, synthetic_frame(angle_deg=angle_deg, seed=seed), pixel_size_nm=pixel_size_nm) \ No newline at end of file + write_synthetic_mrc(path, synthetic_frame(angle_deg=angle_deg, seed=seed, **frame_kwargs), pixel_size_nm=pixel_size_nm) \ No newline at end of file diff --git a/tests/unit/test_frames.py b/tests/unit/test_frames.py index 2a8bafe..1a76f41 100644 --- a/tests/unit/test_frames.py +++ b/tests/unit/test_frames.py @@ -65,3 +65,22 @@ def test_all_buckets_outlier_falls_back_to_overall_consensus(self, tmp_path): assert resolved[paths[0]] == resolved[paths[1]] assert any('no reliable tilt' in str(m) for m in messages) + + + def test_confidence_spike_does_not_skew_overall_consensus(self, tmp_path): + paths, tilt_of = [], {} + # a consistent low-tilt cluster, all striped at 20deg with ordinary confidence + for i, tilt in enumerate([-10, 0, 10]): + path = tmp_path / f'low_{i}.mrc' + write_synthetic_frame(path, angle_deg=20, seed=i) + paths.append(path) + tilt_of[path] = tilt + # a single high-tilt frame with a much sharper peak at a wildly different angle + spike_path = tmp_path / 'spike.mrc' + write_synthetic_frame(spike_path, angle_deg=-60, seed=42, amplitude=600.0, noise_std=1.0) + paths.append(spike_path) + tilt_of[spike_path] = 50 + + resolved = _estimate_per_tilt_angles(paths, tilt_of, angle_outlier_threshold=15.0) + # the low-tilt cluster should still win the overall consensus, not be outvoted by the single spiky frame + assert resolved[tmp_path / 'low_0.mrc'] == pytest.approx(20, abs=1.0) diff --git a/uv.lock b/uv.lock index 10ef717..8664553 100644 --- a/uv.lock +++ b/uv.lock @@ -321,7 +321,7 @@ wheels = [ [[package]] name = "pylisc" -version = "2.2.0" +version = "2.3.0" source = { editable = "." } dependencies = [ { name = "loguru" },