Spread the noiseless signal of a sensor with the kernel of its charge diffusion

[roytsmart]

Stacked on #251 (which is stacked on #244). This is the third part of the diffusion series. It completes the "everything needs to match" goal of #244: the noiseless signal of a sensor now spreads charge with the same kernel that the Monte Carlo samples.

Until now, signal(method="expected") and expose(noise=False) ignored charge diffusion entirely. The Monte Carlo and the uncertainty included it. So the noiseless forward model that inversions use was not the mean of the noisy one.

New

• optika.sensors.kernel_signal(wavelength, *, axis_x, axis_y, ...): the expected electrons in the pixel a photon is absorbed in and in each pixel around it, per absorbed photon. It averages over the photon's position within its pixel and over the depth it is absorbed at, and is a sibling of vmr_signal.
  • It is the quantum yield times the depth average of the differential CCE times the model's per-depth kernel.
  • It sums to the expected electrons per absorbed photon.
• Material:
  • kernel.
  • signal_transposed: the exact adjoint of the noiseless signal.
  • backproject: replaces photons_absorbed.
• Sensor:
  • kernel: the sensor's point-spread function, on axes named after the pixel axes.
  • kernel_variance: the VMR times the kernel. This is the first phase from the plan; an exact variance kernel can come later.
  • expose_transposed: the exact adjoint of expose(noise=False), for gradients.
  • backproject: replaces photons_absorbed.

Changed

• signal(method="expected") with a model of diffusion spreads the photons absorbed in each pixel with kernel_signal.
  • It needs axis_xy and width_pixel, as the Monte Carlo already did.
  • It has the same edge rule, wrap: charge that leaves the grid is lost unless it wraps around.
  • So material.signal(noise=False) and sensor.expose(noise=False) are blurred too.
• The spread is a gather, compiled with numba and parallel over rows. Each output pixel is summed by one thread in a fixed order, so the result is bit-identical for any number of threads (checked on 1, 3, 8 and 48). It allows a different kernel in every pixel.
• The kernel depends on the absorption coefficient, so a direction that varies from pixel to pixel gives every pixel its own kernel. This is the case for collect.
  • Computing a kernel per pixel would cost num² × depth nodes × pixels.
  • Instead, the kernel is computed at 16 optical depths spanning the range and interpolated linearly, one kernel pixel at a time.
  • Against exact per-pixel kernels this agrees to 1e-6 for incidence angles spread over 0 to 30°.
  • With a direction that is uniform over the pixels, as LinearSystem passes it, the kernel is exact.
• backproject is (1/g²) Kᵀ, where g is the gain, the sum of the kernel.
  • Without diffusion it reduces exactly to the old photons_absorbed, e/g.
  • With diffusion, it maps a slowly varying image back onto the photons that produced it. MART and LinearSystem.backproject_from_weights use it.
• AbstractImagingSensor.uncertainty is now private (_uncertainty). Only expose used it.

Breaking

• AbstractImagingSensor.photons_absorbed and AbstractSensorMaterial.photons_absorbed are removed. Use backproject.
• AbstractImagingSensor.uncertainty is removed. Use expose(uncertainty=True).
• AbstractSensorMaterial has three new abstract methods: kernel, signal_transposed and backproject.
• expose(noise=False) changes for any material with a model of diffusion, since its expected image is now blurred.
• signal(method="expected") with a model of diffusion now needs axis_xy and photons that vary along both axes.

Nothing outside optika calls photons_absorbed or uncertainty. ctis calls sensor.expose(noise=False) to probe the gain and assumes the sensor is diagonal, so that needs updating when ctis moves to this release. It is pinned to optika<4.

Things to know

• Gain truncation: the kernel averages over the photons absorbed in the substrate, as the Monte Carlo does, but charge_collection_efficiency assumes an infinitely thick substrate. So for photons the substrate absorbs weakly, the gain with diffusion differs from the gain without it. For the 14 µm CCD97:

  • under 1e-7 from 20 to 5000 Å;
  • 0.06% at 10 Å and 0.05% at 7000 Å;
  • 0.36% at 1 Å and 0.30% at 1 µm.

  The kernel's docstring says so. I left charge_collection_efficiency alone, since its form is deliberate (comparability with Stern 1994).

• Cost: for a 2048 × 1024 frame with a 9 × 9 kernel, per wavelength. The first numpy version of this PR took 0.56 s for one kernel everywhere, and 5.9 s for interpolated per-pixel kernels.

  Kernel	1 thread	8 threads	48 threads
  One kernel everywhere	0.13 s	0.020 s	0.013 s
  Interpolated per pixel	0.33 s	0.053 s	0.032 s

  The transpose costs the same.

• The deferred items from Make models of charge diffusion the one public interface for diffusion #244:

  • cce_backsurface < 1 is now tested against the Monte Carlo.
  • The depth sum needed no chunking: per-pixel directions go through the interpolation, so kernel_signal is evaluated at 16 nodes rather than once per pixel.

Tests

• test_kernel_signal:
  • sum versus the closed-form gain, for a finite and an infinite substrate;
  • with cce_backsurface=1, the result is the quantum yield times kernel_average;
  • without a model, the kernel is a single pixel;
  • num is validated;
  • the errors for missing arguments.
• test_kernel_signal_monte_carlo: the Monte Carlo mean around one pixel equals the kernel, for photons absorbed in the implant with cce_backsurface=0.2, for the Janesick and slab models.
• test_signal_expected_diffusion: a single source spreads with the kernel, with and without wrap.
• test_signal_expected_direction: interpolated per-pixel kernels versus exact ones.
• _blur_test.py: the spread is the product with the explicit matrix whose columns are the kernels, to 1e-14. This covers one kernel everywhere, a kernel per pixel and an interpolated kernel, with and without wrap, forward and transposed. It also checks units, including nodes in a different unit than the positions.
• Material tests:
  • kernel matches the noiseless signal;
  • signal_transposed is the adjoint to 1e-12, including per-pixel directions and wrap;
  • backproject inverts uniform light on a periodic grid to 1e-5.
• Sensor tests, now also on a CCD97 sensor with diffusion:
  • expose of one bright pixel is the kernel;
  • kernel_variance is the VMR times the kernel;
  • expose_transposed is the adjoint of expose;
  • backproject inverts expose without diffusion, and inverts uniform light in the middle of a wide sensor with it.
• Locally: the full suite passes (23,923 tests), and the patch is fully covered (436 of 436 lines).

🤖 Generated with Claude Code

https://claude.ai/code/session_01BAw6Lv7Cy2ULGmsMJ4b1dE

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (5ac92e7) to head (60a2c44).

Additional details and impacted files

@@                    Coverage Diff                    @@
##           feature/diffusion-slab      #255    +/-   ##
=========================================================
  Coverage                  100.00%   100.00%            
=========================================================
  Files                         131       133     +2     
  Lines                       11074     11553   +479     
=========================================================
+ Hits                        11074     11553   +479     

Flag	Coverage Δ
unittests	100.00% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[roytsmart]

Request: the covariance between pixels from charge diffusion (motivated by Hinode/XRT)

This stack covers almost everything needed to model the noise of the Hinode/XRT CCD. The missing piece is the covariance between neighboring pixels. I'd suggest it as a follow-up PR stacked on this one rather than more scope here.

Motivation

hinode (sun-data/hinode v0.2.0) computes the per-pixel uncertainty of XRT as compound-Poisson photon noise, $\sigma^2 = F,S$ with $F = \langle s^2\rangle / \langle s\rangle$ (in DN per photon), plus the Kobelski et al. (2014) uncertainty map. Like Kobelski, it assumes each photon's charge lands in one pixel.

I measured the real noise from differences of consecutive frames:

• Data: XRT Al_poly from the ESIS flight, hinode.xrt.open(time_start="2019-09-30T18:06:11", time_stop="2019-09-30T18:11:01", filter="Al_poly", leak=True, uncertainty=True). That is 10 pairs of consecutive 16.385 s frames, 18.5 s apart, JPEG quality 90, full resolution, about 1.26 million pixel differences.

• Pixels left out: the jet, saturated pixels, bright features, and pixels where the 0.02–0.27″ pointing change would move structure. The model signal has the visible-light leak subtracted.

• Variance: the measured variance is 0.15–0.19 of the model, depending on the assumed F (2.98 or 2.47 DN/ph, for log T 6.3 or 6.2). It rises slowly with signal. So σ is overestimated about 2.4×.

• Correlation of the noise between neighbors:

  pair of pixels	correlation
  side by side, inside an 8 × 8 JPEG block	+0.41 to +0.45
  side by side, across a block edge	+0.27
  diagonal	+0.13
  two pixels apart	+0.05

Interpretation. The across-edge, diagonal and two-apart values match a charge cloud of about 5 µm on XRT's 13.5 µm pixels, as in the ccd-snr-paper. That spread puts two electrons from one photon in the same pixel only about 36% of the time, which alone would cut the variance to about 0.36. The extra correlation inside blocks, and the ratio rising with signal, are JPEG quantization zeroing the high-frequency DCT components of the noise. That accounts for roughly another ×0.45.

I couldn't find this in the XRT literature. Kano et al. (2008) give the CCD (e2v CCD42-40, back-illuminated, 13.5 µm, 57 e⁻/DN) but no depletion depth and nothing on charge spreading. Kobelski et al. (2014) estimate the photon noise from a model and never compare it with measured scatter.

What it changes. In a joint AIA + XRT DEM (utu.dem.plowman), the overestimated σ lets the inversion stop at reduced χ² = 1 with XRT at +1.75σ. With the measured noise, the emission measure above log T 6.6 drops about 37% in the quiet Sun.

Plan. Model the XRT CCD with this stack, simulate XRT's JPEG in hinode, and run the same frame-difference analysis on simulated frames. Then compare the variance ratio and the correlation pattern with the data. The correlations are the cleanest handle on thickness_depletion, cleaner than the variance, which JPEG also lowers.

What this stack already covers

• The CCD material: BackIlluminatedSiliconSensorMaterial.fit_eqe, with Kano et al. (2008)'s QE of XRT: 0.93 at 13 Å, 0.61 at 45 Å, 0.46 at 114 Å, 0.56 at 304 Å.
• Diffusion: JanesickDiffusionModel or SlabDiffusionModel, with thickness_depletion as the free parameter.
• Per-pixel variance with diffusion: vmr_signal, through probability_same_pixel, for X-rays from 1 Å up.
• The mean point-spread kernel: kernel_signal.
• Simulated flat fields: signal(method="monte-carlo", diffusion=..., axis_xy=..., wrap=True), to feed a JPEG simulation.
• Al_poly's spectrum: the caller evaluates along a wavelength axis and weights by the electrons each wavelength contributes. No change is needed here for that.

What's missing: the two-point covariance

For uniform illumination, the photon arrivals are compound Poisson. So the covariance between pixels offset by $\Delta$ is the absorbed-photon rate times $E[X_0 X_\Delta]$ for a single photon, the same faint-exposure argument as the paper. A photon absorbed at depth $z$ delivers $m$ surviving electrons. They are spread over pixels with probabilities $q_j$, with $\sum_j q_j = 1$, so

$$E[X_j X_k] = E[m],q_j,\delta_{jk} + E[m(m-1)],q_j,q_k .$$

Summing over where the photon lands in the uniform field, and averaging over its position within the pixel and over depth, gives the covariance per absorbed photon, in electrons²:

$$\mathrm{Cov}(\Delta) = \langle n\rangle \langle\eta\rangle, \delta_{\Delta 0} + \langle n\rangle \left(\langle n\rangle + \mathcal{F} - 1\right) \left\langle \mathcal{P}_\Delta, \eta^2 \right\rangle, \qquad \mathcal{P}_\Delta(z) = \Big\langle \sum_j q_j, q_{j+\Delta} \Big\rangle_{\text{position in pixel}} .$$

Here $\mathcal{P}_0 = $ probability_same_pixel, and the outer $\langle\cdot\rangle$ is the absorption-weighted depth average of average_depth. At $\Delta = 0$ this is exactly what vmr_signal computes, times the mean electrons per absorbed photon.

kernel_variance (VMR × kernel) spreads variance, but it isn't this covariance. Your PR description already flags an exact variance kernel as later work, and this would be it.

Proposed API

1. AbstractDiffusionModel.kernel_pair(depth, thickness_substrate, width_pixel, axis_x, axis_y, num=None)
   • It returns $\mathcal{P}_\Delta(z)$ on the same offset grid as kernel: the probability that two electrons created at depth by the same photon are collected in pixels offset by $(\Delta x, \Delta y)$, averaged over the photon's position within its pixel.
   • For a Gaussian cloud it separates by axis: $\mathcal{P}\Delta = \mathcal{P}^x{\Delta x},\mathcal{P}^y_{\Delta y}$, with $\mathcal{P}^x_\Delta = \int_0^1 du \sum_k q_k(u), q_{k+\Delta}(u)$.
   • For SlabDiffusionModel each electron draws its own mixture component, so it sums over pairs of components exactly as probability_same_pixel already does.
   • Its center must equal probability_same_pixel.
2. optika.sensors.covariance_signal(wavelength, *, axis_x, axis_y, direction, n, n_substrate, thickness_implant, thickness_substrate, diffusion, width_pixel, cce_backsurface, temperature, num=None)
   • It returns $\mathrm{Cov}(\Delta)$ per absorbed photon as a FunctionArray over offsets, a sibling of vmr_signal and kernel_signal.
   • The $\eta^2$ weighting has to sit inside the depth average, as in vmr_signal. So it belongs at the material level, not in the model.
3. Material and sensor methods, e.g. AbstractSensorMaterial.kernel_covariance and AbstractImagingSensor.kernel_covariance, mirroring kernel and kernel_variance.

Tests that would pin it down

• Center against vmr_signal: the center of covariance_signal equals vmr_signal × the mean electrons per absorbed photon.
• Sum against no diffusion: the sum over all offsets equals the variance without diffusion, vmr_signal(diffusion=None) × the same mean, within the kernel's truncation. Diffusion moves charge between pixels without creating or destroying it.
• No model: it is a single pixel equal to the undiffused variance.
• kernel_pair against probability_same_pixel: its center equals probability_same_pixel, and its sum is 1 within truncation.
• Monte Carlo: the empirical covariance at offsets 0, (1, 0), (1, 1) and (2, 0) of flat fields from signal(method="monte-carlo", diffusion=..., wrap=True) matches, for the Janesick and slab models, including cce_backsurface < 1.

Not for optika

The JPEG simulation is XRT-specific and will live in hinode, with XRT's quantization tables per QTABLE.

[roytsmart]

Re #255 (comment): implemented in draft #256, stacked on this PR.

• The API, as proposed:
  • AbstractDiffusionModel.kernel_pair is built from the new public mixture (Add SlabDiffusionModel, a model of diffusion across the field-free region #251). Each pair of Gaussians gives _kernel_1d with variance vᵢ + vⱼ.
  • optika.sensors.covariance_signal gives the covariance in e²/ph per absorbed photon.
  • kernel_covariance is on the materials and on AbstractImagingSensor.
• One change to the formula: covariance_signal is the undiffused variance at the center plus a term that only moves variance between pixels. So its center equals vmr_signal × the gain and its sum equals the undiffused variance exactly, at every wavelength from 1 to 7000 Å. The direct form would mix vmr_signal's infinite substrate with the truncated depth average.
• Tests: all the tests you listed are in.
  • The Monte Carlo test checks offsets 0, (1, 0), (1, 1) and (2, 0) for both models with cce_backsurface = 0.2. Its tolerance is six times the standard error from the scatter between 64 frames.
  • For a CCD97 with 16 µm pixels at 304 Å, the correlation is 0.239 side by side and 0.084 diagonally.
• Caveat for fitting XRT: the model, and the Monte Carlo, treat each electron as diffusing independently from one point. At X-ray energies the initial size and self-repulsion of the charge cloud correlate the electrons and narrow the pair separations. A fitted thickness_depletion would absorb that.

🤖 Generated with Claude Code