mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 80 Wiki Activity

Scaling: a correction surface is fitted on the negative observations too

Browse Source 
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 <I>/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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
leonarski_f
2026-09-21 16:14:30 +02:00
co-authored by Claude Opus 5
parent 035a58d9e0
commit 1d1a40d908
1 changed files with 32 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+32 -15
View File
@@ -3174,19 +3174,20 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &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<int32_t> &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<int32_t> &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<double>(o.I) * o.corr * a, sc = static_cast<double>(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<int32_t> &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<double> 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;