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;