From 1d1a40d908c935b1a85edcb4bfd23ef4bc98bab2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 21 Sep 2026 16:14:30 +0200 Subject: [PATCH] Scaling: a correction surface is fitted on the negative observations too The per-cell fit in ApplyCellSurface dropped every observation with Is <= 0. Both of its sums are linear in Is - there is no logarithm and no division by it - so the filter protected nothing arithmetic; the only hazard of keeping the negatives is a cell whose `cross` sums to <= 0, which now keeps its factor for that round instead of dividing by it. What the filter did was select on the noise: it kept the up-fluctuated half of a weak cell, inflating `cross` by the truncated mean, so each cell was pulled down by an amount set by its own signal-to-noise. Per round, on a Wilson-distributed cell with six equivalents: 0.1% at /sigma = 5, 1.0% at 2, 4.9% at 1, 15% at 0.5. That push is almost entirely along the surface's resolution null direction, which is why the surfaces walked. Simulated with a flat truth and a 12% cos(2 phi) flat field: - no gauge: shell at I/sigma = 1 suppressed 40% after 30 rounds with the filter, flat to 1% without it; - 16x16 detector grid, gauge on, outer I/sigma = 0.5: with the filter the fitted field carries 1.9x its true amplitude at CC 0.48 to the truth and the gauge removes a 5x ramp; without it 1.04x, CC 0.97, 1.09x. On 12 open rotation sets against the gauge build (rc172-2): - the null walk the gauge removes collapses where it was large (45.9x -> 2.0x, 5.0x -> 2.4x, 3.1x -> 1.7x; a 56x and a 4.9e5x walk are gone because those surfaces no longer pass cross-validation); - held-out gains lose the part the filter inflated (23.2 -> 20.1%, 14.1 -> 12.9%); - sets at I/sigma >= 10: per-shell CC(Fo, Fc) against the deposited model moves by -0.006..+0.003 (fine half), ISa unchanged; - the two weakest sets (I/sigma 6.8 and 4.6) lose a surface the gate no longer admits, and with it 0.019 / 0.031 of fine-half CC(Fo, Fc). Forcing the unfiltered surface on the first recovers all of it, so there it is the gate that refuses a real correction; on the second the unfiltered surface is worth nothing and the filtered one's benefit is not a correction the fit can reproduce. "NOT settled in 30 rounds" does not go away (10 of 17 applied surfaces -> 8 of 15): it has a second cause the filter was not driving. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/RotationScaleMerge.cpp | 47 +++++++++++++------ 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index a56ee38dc..67a8582b8 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3174,19 +3174,20 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // the same on every round, so the surface walks along it at constant speed until the [0.25, 4] // clamp stops it. That is what "NOT settled in n round(s)" reports, and why raising the cap from // 3 to 30 rounds took the walk ten times further rather than converging it. - // Something has to push it: here the `Is > 0` filter in the fit below, which drops the negative - // half of a weak cell's observations and leaves the survivors up-fluctuated, so `cross` is - // inflated where the signal is weak - a bias that is a function of the cell's signal-to-noise, - // hence of resolution. Simulated on synthetic data with a flat truth: an outer shell at - // I/sigma = 1 is suppressed 4% after 3 rounds and 24% after 30, while the even/odd held-out gain - // that gates the surface reads 0.1% - the cross-validation cannot see it either, because both - // halves carry it equally. - // So it is removed at the parameterisation instead of being bounded by a round count: after each - // round the surface is divided by its own weighted geometric mean WITHIN each resolution shell, - // which projects out the component that is a function of d and leaves the azimuthal / positional - // variation - what absorption and a flat field physically are. On the same simulation that - // recovers a 12% cos(2 phi) flat field to CC 0.986 at full amplitude while the radial ramp goes - // to 1.00. The projection removes the overall level too, so it subsumes the global gauge below. + // Something has to push it, and what used to push it was a `Is > 0` filter on the observations + // entering the fit: it dropped the negative half of a weak cell's measurements and left the + // survivors up-fluctuated, so `cross` was inflated where the signal is weak - a bias that is a + // function of the cell's signal-to-noise, hence of resolution. That filter is gone (see the fit + // below), which takes the bulk of the push away at the source; the gauge stays because the null + // direction is exact, so whatever else pushes along it - finite-sample noise in the update, the + // clamp - would still walk it. + // The null direction is removed at the parameterisation rather than being bounded by a round + // count: after each round the surface is divided by its own weighted geometric mean WITHIN each + // resolution shell, which projects out the component that is a function of d and leaves the + // azimuthal / positional variation - what absorption and a flat field physically are. On the + // same simulation that recovers a 12% cos(2 phi) flat field to CC 0.986 at full amplitude + // while the radial ramp goes to 1.00. The projection removes the overall level too, so it + // subsumes the global gauge below. // Equal-occupancy shells in s^2 = 1/d^2. The count is second order (8, 20 and 40 shells all leave // the residual ramp inside 1.5%); 20 is fine enough to follow a curved drift and coarse enough // that each shell still holds tens of thousands of observations. @@ -3300,6 +3301,15 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // the gauge fix spreads over the whole surface and the n_iter rounds compound, and which no // even/odd cross-validation can see because both halves carry it equally. The per-frame // scale (FitPerFrameG) already regresses this way round. + // A NEGATIVE observation is a legitimate measurement of a weak reflection and enters both + // sums as it stands. Both are linear in Is - there is no logarithm and no division by it + // here, unlike the relative-B fits above, which regress log(Iref/Is) and do need Is > 0 - + // so the only thing a sign filter would protect is the division by `cross` below, and that + // is cheaper to guard directly. Selecting on the sign instead keeps the up-fluctuated half + // of a weak cell and inflates `cross` by the truncated mean, which is a bias in the cell's + // own signal-to-noise: simulated on a flat truth, it suppressed a shell at I/sigma = 1 by + // 40% over 30 rounds and left a 16x16 flat field at outer I/sigma = 0.5 with 1.9x its true + // amplitude and CC 0.48 to the truth, where keeping the negatives gives 1.04x and 0.97. // Chunked over the subset, each thread summing into its own cells: unlike the pass above // there is no ordering that keeps threads off each other's bins, and there are only ncell // of them, so per-thread copies are cheap and the fixed chunking keeps it reproducible. @@ -3314,7 +3324,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int if (sw[o.group] <= 0.0) continue; const double Iref = swI[o.group] / sw[o.group], a = A[o.cell]; const double Is = static_cast(o.I) * o.corr * a, sc = static_cast(o.sigma) * o.corr * a; - if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0) || !(sc > 0.0)) continue; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(sc > 0.0)) continue; const double w = 1.0 / (sc * sc); xcross[o.cell] += w * Is * Iref; xref2[o.cell] += w * Iref * Iref; } @@ -3326,7 +3336,14 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int std::nth_element(dsorted.begin(), dsorted.begin() + dsorted.size() / 2, dsorted.end()); const double lambda = 0.1 * std::max(1e-30, dsorted[dsorted.size() / 2]); std::vector upd(ncell, 1.0); - for (int c = 0; c < ncell; ++c) upd[c] = (ref2[c] + lambda) / (cross[c] + lambda); + // A cell whose observations carry no net signal can sum to `cross` <= 0 now that the + // negatives are in, and the Tikhonov pull does not always lift it clear of zero. Such a + // cell has nothing to say about its own factor, so it keeps the one it has rather than + // dividing by a number that can be near zero or negative - which would reach the log below + // as a NaN or as a factor pinned on the clamp and, through the shell means, poison every + // other cell. It is the same `cross > 0` the gauge and the step size already test. + for (int c = 0; c < ncell; ++c) + if (cross[c] > 0.0) upd[c] = (ref2[c] + lambda) / (cross[c] + lambda); double logsum = 0.0, wsum = 0.0; for (int c = 0; c < ncell; ++c) if (cross[c] > 0.0) { logsum += cross[c] * std::log(upd[c]); wsum += cross[c]; } const double gm = wsum > 0.0 ? std::exp(logsum / wsum) : 1.0;