FF files: store the average at full 16-bit precision; fix average rounding and stdpixel variance

[Cybis320]

Motivation

The FF average plane discards most of the information the 256-frame mean actually contains, and the discarded fraction corrupts more than just the average:

• The trimmed mean of 252 frames has intrinsic noise of stdpixel/sqrt(252) ≈ 0.2–0.4 ADU on real stations (measured on USV001 and ES001D archives) — at or below the 1 ADU quantization step. On USV001, the majority of sky pixels are quantization-dominated.
• The stored average is floor-truncated (C integer division), biasing it −0.5 ADU. Meteor photometry sums of maxpixel − avepixel run correspondingly hot.
• The quantization error is spatially coherent (flat sky floors to the same integer everywhere), so background estimators can't remove it: measured faint-star aperture photometry scatter is ×1.7–×5.2 worse than the information limit, i.e. 0.55–1.8 mag of point-source depth on the average image is lost. The sub-ADU information demonstrably survives 8-bit frame quantization (dithered by per-frame noise); it is destroyed only by the output cast.
• The floored mean is reused in the variance term (acc*mean instead of acc²/n), inflating stdpixel by an extra ≈ mean × frac(mean) ADU² per pixel: measured +39–48% at typical sky levels and up to +288% on quiet cameras. The old NOTE claiming the non-standard calculation "does not matter" was quantitatively false.

Changes

1. Full-precision average plane. compressFrames additionally returns the trimmed mean in 8.8 fixed point (uint16, units of 1/256 ADU, rounded). With ff_avepixel16: true (new config option, FITS only) it is written as the AVEPIXEL HDU with an AVEFRAC=8 header keyword. FFfits.read keeps it in ff.avepixel16 and derives the legacy 8-bit ff.avepixel view by rounding off the fractional bits, so every existing consumer keeps working unchanged. Legacy files, .bin, and native 16-bit camera files are unaffected (covered by tests).
2. Rounded 8-bit average. The uint8 plane is now (ave16 + 128) >> 8 — zero-mean rounding (measured bias +0.01 ADU vs −0.494 before), one shared derivation between compressor and readers.
3. Correct stdpixel. The variance uses the proper sample formula computed in double precision, and the standard deviation is rounded. Residual error vs the true trimmed std is only the uint8 output rounding (−1…−12%, worst for quiet cameras where σ≈1.5).
4. Precision consumers. extractStarsFF (CALSTARS) and SkyFit2's three measurement paths (star extraction, centroiding, manual-photometry background) prefer avepixel16 when present, with identical fallback otherwise.
5. Incidental fix: modern astropy raises on explicitly memory-mapped BZERO-scaled HDUs, which silently broke reading any native 16-bit FF file. FFfits.read now probes and falls back to a plain read for scaled files (all planes are copied out regardless).

Measured results

• End-to-end through the production write/read path: faint-star aperture scatter improves ×5.1–5.6; +1.8 mag of faint-end depth on the average image at SNR=5.
• SkyFit2 photometry background median: −0.23 ADU error on the 8-bit path vs 0.00 on the full-precision path (synthetic scene).
• Compressor outputs validated bit-exact against a numpy reference, including saturation, ties, and low-noise edge cases (Tests/TestFFAvepixel16.py).
• On files without the new plane, star extraction output is bit-identical to the previous code (verified field-for-field on real archived FFs), and FrameInterface, thresholdFF, MakeFlat, StackFFs, reconstructFrame, and adjustLevels produce identical results between new-format and legacy copies of the same data.

Cost and compatibility

• FF size: 3.70 → 4.62 MB at 720p (+25%). Composes with the RICE compression of FITS Handling: Introduce Lossless Compression for HDUs #278: RICE + full-precision average measures 2.63 MB — still ~29% smaller than today's uncompressed files.
• Mixed per-HDU dtypes are core FITS; a converted real FF (BITPIX 8/8/16/8) round-trips losslessly through astropy and reads cleanly with cfitsio.
• Behavioral changes to test (why this is a draft):
  • The rounded 8-bit average shifts values by 0/+1 ADU relative to the old floor.
  • Corrected stdpixel drops ~30–40% from the old inflated values, so k1-based detection thresholds become effectively more sensitive; field testing will show whether k1/j1 need a compensating migration or the correct values can stand.
  • Readers older than this branch will misread the uint16 average plane on new files; ff_avepixel16: false restores fully legacy output.

Update: correctness of the average/std computations themselves, and an audit of the same error class

Beyond storing more bits, the computations had four systematic errors, all fixed here (each validated against exact numpy references and, where noted, measured):

1. Floor instead of round (C division semantics): average biased −0.5 ADU. Now rounded (bias +0.01).
2. stdpixel variance formula: the truncated mean was reused in the acc*mean term instead of acc²/n, inflating stdpixel by ~mean*frac(mean) ADU² — +39–48% at typical sky, up to +288% on quiet cameras. Now the correct sample variance in double precision. Note: corrected stdpixel is ~30–40% lower, so k1-based thresholds become effectively more sensitive — the main field-testing question of this draft.
3. Unbalanced trim: removing only the top 4 frames biases the mean −0.04 σ (−0.18 ADU at σ=4.5). The bottom 4 are now trimmed as well (measured bias +0.001 ADU).
4. Averaging in the encoded domain: the mean of gamma-encoded samples is biased low (Jensen; measured −0.24 ADU at γ=0.6, sky 40). With a configured camera gamma ≠ 1, frames are decoded via LUT, averaged in the linear domain, and re-encoded — the plane stays gamma-encoded so all consumers are unaffected, and their own gamma correction now recovers the true linear mean. AVEGAMMA records the gamma used; γ=1 keeps the exact integer path.

A codebase audit for the same class (floor-instead-of-round, premature narrow-dtype truncation) found and fixed:

• FFMimickInterface (video/image-sequence inputs) replicated both compressor bugs: −0.48 ADU average bias, stdpixel +26–146% → detection/photometry on non-FF inputs. Fixed (measured −0.02 ADU, −1%), and subsequently brought to full parity with the compressor: symmetric max/min trim, linear-domain averaging with the configured gamma (measured +0.007 ADU vs −0.202 encoded-domain through the production FrameInterface path), and an avepixel16 fixed-point plane for 8-bit content, so the SkyFit2 precision paths work on video inputs with no further changes. The gamma white point and fixed-point gate come from config.bit_depth (known at construction) rather than the late-bound output dtype, and the gate checks the data range so sum-binned frames cannot wrap.
• applyFlat and gammaCorrectionImage floored on the cast back to integer types (−0.5 ADU per application).
• binImage: floored 'avg'; float flats truncated to uint16 before binning; uint8 'sum' binning wrapped mod 256.
• MakeFlat: flats floor-truncated (a field-dependent multiplicative photometric error, since the flat divides every image); >8-bit darks silently wrapped modulo 256.
• Truncated report-time casts: CALSTARS intensity/amplitude/background, FTPdetectinfo level/background, detection & SkyFit2 intensity sums; Astra float32→int32 background subtraction; FR-file avepixel reconstruction.

Reported but deliberately not changed here (each needs its own validation): preprocessFF dark/flat-corrects maxpixel and avepixel but never stdpixel (both FF and non-FF paths), so k1*stdpixel thresholds are mis-scaled wherever the flat deviates from unity (~30% too sensitive in strongly vignetted corners); applyDark's cv2.subtract clips negatives for integer inputs only, so meteor sums currently depend on whether gamma correction is active; FFMimickInterface still averages in the encoded domain. (An earlier revision of this description claimed non-FF detection subtracts an uncalibrated average from calibrated frames - that was wrong: preprocessFF runs after every chunk load.)

🤖 Generated with Claude Code