GridBased/Notebooks/KnifeEdge_2filters.ipynb

Knife edge, 2 filters — GridBased against the published chains¶

Same measurement, same ROI, same ESF construction and same erf fit as DeepLearning/Notebooks/KnifeEdge_2filters.ipynb, so the only thing that differs between the rows of the final table is the reconstruction itself. The flat field sits at 1.81 % occupancy inside the ESF ROI, 2.4x the 3-filter point.

Edge width (revised 2026-09-28). The ESF is binned and fitted exactly as Plotting/fig_knife_edge.py --fit odr does it: generate_presampled_esf() and get_odr_fitting() from DeepLearning/src/funcs_postprocess.py (bin content Σk/Σf), symmetric ±4σ window, ESF bin width 0.3 µm. The ODR model is p0/2 (1 + erf((x - p1)/(sqrt(2) sigma))) + p3, so the fit gives sigma, the width of the Gaussian LSF, which is what the tables quote, with MTF50 [lp/mm] = 1000 sqrt(ln 2 / 2) / (pi sigma); MTF50 and MTF(200 lp/mm) are also measured from the LSF (Plotting/knife_edge_numbers.py: first crossing, central-difference correction), and the relative DQE uses the latter. Only the flattened GridBased reconstruction is analysed. Same processing as KnifeEdge_poly2to20keV_Moench040.ipynb.

Occupancy is that of the flat field inside the ESF ROI: the reconstructed photons per 25 µm pixel per frame over ROI // 10 (x 33–46, y 26–49 of the 99 × 99 frame), the pixels over which DeepLearning/Notebooks/KnifeEdge_3filters_pos0.ipynb sums the deposited energy (0.705 % at 3 filters from the energy, 0.747 % here from the photons of the 12 keV model).

Until 2026-09-28 this notebook took the fitted sigma for sqrt(2) sigma (so its MTF50 was sqrt(2) too high), used 0.5 µm bins, and binned the ESF with an older copy of generate_presampled_esf (inverse-variance mean of k/f, biased low in the toe, sigma ~3 % too small). The previous version is in Backup_260928_pre_sigma_fix/.

The ROI, the oversampled ESF and the erf fit are the ones from DeepLearning/Notebooks/KnifeEdge_3filters_pos0.ipynb; the edge sits at the same place at all three flux points, so the same ROI holds throughout and the three notebooks are directly comparable.

Run in the pytorch_nightly environment (ROOT + scipy + jupyter).

In [1]:

import json
import sys
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

REPO = Path('/home/xie_x1/MLXID/GridBased')
sys.path.insert(0, '/home/xie_x1/MLXID/DeepLearning/src')   # the code Plotting/ uses
sys.path.insert(0, '/home/xie_x1/MLXID/Plotting')
from ROOT import TCanvas, gROOT
gROOT.SetBatch(True)
from funcs_postprocess import generate_presampled_esf, get_odr_fitting
from knife_edge_numbers import mtf as mtf_from_lsf           # first crossing, sinc-corrected

GRID = Path('/mnt/sls_det_storage/moench_data/MLXID/InferenceResults/GridBased')
PUBLISHED = Path('/home/xie_x1/MLXID/DeepLearning/InferenceResults')
ETA = Path('/home/xie_x1/MLXID/EtaInterpolation/Results')

EDGE_RUN, FLAT_RUN = 'Edge_2filters_12keV', 'Flat_2filters_12keV'
BIN_UM = 2.5                              # super-pixel = 25 um / 10
ROI = (330, 465, 260, 490)                # x1, x2, y1, y2, as published
CROP = slice(10, -10)                     # 1010 -> 990, as published
BIN_WIDTH_UM = 0.3                        # presampled ESF bin, as Plotting/ and the manuscript

plt.rcParams.update({'figure.dpi': 110, 'image.cmap': 'gray',
                     'image.origin': 'lower', 'image.interpolation': 'nearest'})

Welcome to JupyROOT 6.28/10

1. Load every reconstruction of the same measurement¶

In [2]:

def crop(image):
    """The published frames are cropped by 10 super-pixels on each side."""
    return image[CROP, CROP] if image.shape[0] == 1010 else image


for run in (EDGE_RUN, FLAT_RUN):
    if not (GRID / run / 'summary.json').exists():
        raise SystemExit(
            f'{run} has not been reconstructed yet.  Run:\n'
            f'  python scripts/05_infer_measurement.py '
            f'--config Configs/infer_measurement_12keV.yaml --run {run}')

grid_summary = json.load(open(GRID / EDGE_RUN / 'summary.json'))
flat_summary = json.load(open(GRID / FLAT_RUN / 'summary.json'))


def grid_image(run, flattened):
    name = 'binned_image_flattened.npy' if flattened else 'binned_image.npy'
    return np.load(GRID / run / name).astype(float)


variants = {
    'eta':            (np.load(ETA / '2603MaxIV_Edge2Filters_12keV_fromClusters.npy'),
                       np.load(ETA / '2603MaxIV_Flat2Filters_12keV_fromClusters.npy')),
    'GridBased flat': (grid_image(EDGE_RUN, True), grid_image(FLAT_RUN, True)),
}

EXPERIMENT = {1: '1ph_260512_12keV_v260511_06',
              2: '2ph_260610_12keV_v260610_01',
              3: '3ph_260611_12keV_v260611_00'}
KNIFE, FLAT = 'KnifeEdge_2filters', 'FlatField_2filters'


def published_frame(run, n_photons):
    name = f'{n_photons}Photon_ML_superFrame_etaInterpolated.npy'
    return np.load(PUBLISHED / run / EXPERIMENT[n_photons] / 'augX1' / name).astype(float)


cluster_edge = {n: published_frame(KNIFE, n) for n in (1, 2, 3)}
cluster_flat = {n: published_frame(FLAT, n) for n in (1, 2, 3)}
variants.update({
    'cluster 1ph':   (cluster_edge[1], cluster_flat[1]),
    'cluster <=2ph': (cluster_edge[1] + cluster_edge[2], cluster_flat[1] + cluster_flat[2]),
    'cluster <=3ph': (cluster_edge[1] + cluster_edge[2] + cluster_edge[3],
                      cluster_flat[1] + cluster_flat[2] + cluster_flat[3]),
})

ORDER = ['eta', 'cluster 1ph', 'cluster <=2ph', 'cluster <=3ph', 'GridBased flat']
variants = {name: (crop(e), crop(f)) for name, (e, f) in variants.items()}
variants = {name: variants[name] for name in ORDER}

def flat_occupancy(run):
    # reconstructed photons per 25 um pixel per frame in the flat field, inside the
    # ESF ROI: ROI // 10 in pixels of the 99 x 99 frame, both ends included as in
    # DeepLearning/Notebooks; +1 because the 101 x 101 frame is cropped by 1 px
    summary = json.load(open(GRID / run / 'summary.json'))
    native = np.load(GRID / run / 'native_image.npy') / summary['n_frames']
    x0, x1, y0, y1 = (v // 10 for v in ROI)
    return float(native[y0 + 1:y1 + 2, x0 + 1:x1 + 2].mean())


occupancy = flat_occupancy(FLAT_RUN)
print(f"GridBased model epoch {grid_summary['epoch']}, "
      f"{grid_summary['n_frames']:,} frames, "
      f"occupancy {occupancy*100:.3f} % (flat field, ESF ROI)")
print(f"  sub-pixel excess rms {grid_summary['subpixel_nonuniformity_before']*100:.1f} % -> "
      f"{grid_summary['subpixel_nonuniformity_after']*100:.1f} % after flattening\n")
print(f"{'variant':16s} {'edge photons':>14} {'flat photons':>14}  relative")
reference = variants['cluster <=3ph'][0].sum()
for name, (e, f) in variants.items():
    print(f'{name:16s} {e.sum():14.4e} {f.sum():14.4e}  {e.sum()/reference:8.3f}')

GridBased model epoch 135, 7,000,000 frames, occupancy 1.808 % (flat field, ESF ROI)
  sub-pixel excess rms 21.6 % -> 3.8 % after flattening

variant            edge photons   flat photons  relative
eta                  1.5670e+08     2.7062e+08     0.710
cluster 1ph          1.5670e+08     2.7062e+08     0.710
cluster <=2ph        2.0504e+08     3.5599e+08     0.929
cluster <=3ph        2.2076e+08     3.8444e+08     1.000
GridBased flat       2.3053e+08     4.0162e+08     1.044

In [3]:

fig, axes = plt.subplots(1, len(ORDER) if len(ORDER) < 4 else 3, figsize=(15, 4.6))
extent = [0, 990 * BIN_UM, 0, 990 * BIN_UM]
shown = ORDER if len(ORDER) < 4 else ['eta', 'cluster <=3ph', 'GridBased flat']
for ax, name in zip(np.atleast_1d(axes), shown):
    image = variants[name][0]
    im = ax.imshow(image, extent=extent, vmax=np.percentile(image, 99.5))
    ax.add_patch(plt.Rectangle((ROI[0]*BIN_UM, ROI[2]*BIN_UM),
                               (ROI[1]-ROI[0])*BIN_UM, (ROI[3]-ROI[2])*BIN_UM,
                               edgecolor='red', facecolor='none', lw=1.4))
    ax.set_title(name); ax.set_xlabel('x [um]'); ax.set_ylabel('y [um]')
    plt.colorbar(im, ax=ax, fraction=0.046)
fig.suptitle('knife edge, with the published ESF ROI')
fig.tight_layout()

2. Edge spread function¶

generate_presampled_esf fits an erf to every row of the ROI, fits a straight line through the edge positions, and re-samples the ROI onto a common axis; get_odr_fitting then fits the erf with x and y errors inside |x - p1| <= 4 sigma (the fit of Plotting/fig_knife_edge.py --fit odr).

In [4]:

def analyse(edge_image, flat_image, label):
    esf = generate_presampled_esf(edge_image, flat_image, ROI,
                                  pixel_size_um=BIN_UM, bin_width_um=BIN_WIDTH_UM,
                                  chi2_threshold=15.0, residual_tolerance_px=1.0)
    x, z = esf['numpy']['x'], esf['numpy']['z']
    odr = get_odr_fitting(x, z, esf['numpy']['x_err'], esf['numpy']['z_err'], n_sigma=4.0)
    sigma, sigma_err = odr['sigma_um'], odr['sigma_err_um']     # sigma of the Gaussian LSF
    try:
        mtf50_lsf, mtf200_lsf = mtf_from_lsf(dict(x=x, z=z, pars=odr['beta']))
    except Exception:
        mtf50_lsf = mtf200_lsf = np.nan
    n_fit = int(odr['mask'].sum())
    ndf = max(n_fit - len(odr['beta']), 1)
    x_fit = x[odr['mask']]
    return dict(label=label, x=x, z=z, sigma=sigma, sigma_err=sigma_err, sigma_gauss=sigma,
                chi2_ndf=odr['chi2_ndf'], chi2=odr['chi2_ndf'] * ndf, ndf=ndf,
                mtf50=1000 * np.sqrt(np.log(2) / 2) / (np.pi * sigma),
                mtf10=1000 * np.sqrt(np.log(10) / 2) / (np.pi * sigma),
                mtf50_lsf=mtf50_lsf, mtf200_lsf=mtf200_lsf,
                rows_kept=esf['valid_lines_ratio'],
                edge_slope=esf['edge_slope'], edge_centre=odr['beta'][1],
                fit_points=n_fit, fit_range=(x_fit.min(), x_fit.max()),
                photons=float(edge_image.sum()), odr=odr)


results = {name: analyse(edge, flat, name) for name, (edge, flat) in variants.items()}

print(f"{'variant':16s} {'sigma [um]':>16} {'chi2/ndf':>9} {'MTF50':>7} {'MTF50 LSF':>9} "
      f"{'MTF(200) LSF':>12} {'rows':>6}")
for name, r in results.items():
    print(f"{name:16s} {r['sigma']:7.3f} +- {r['sigma_err']:.3f} {r['chi2_ndf']:9.2f} {r['mtf50']:7.1f} "
          f"{r['mtf50_lsf']:9.1f} {r['mtf200_lsf']:12.3f} {r['rows_kept']*100:5.0f}%")

Final Edge Midpoint Fit: slope = -0.0858, intercept = 432.3990, chi2/ndf = 0.0751
Final Edge Midpoint Fit: slope = -0.0854, intercept = 432.2740, chi2/ndf = 0.0576
Final Edge Midpoint Fit: slope = -0.0855, intercept = 432.4402, chi2/ndf = 0.0479
Final Edge Midpoint Fit: slope = -0.0850, intercept = 432.2537, chi2/ndf = 0.0406
Final Edge Midpoint Fit: slope = -0.0858, intercept = 432.5887, chi2/ndf = 0.0361
variant                sigma [um]  chi2/ndf   MTF50 MTF50 LSF MTF(200) LSF   rows
eta                1.483 +- 0.024      1.61   126.4     131.7        0.187   100%
cluster 1ph        1.352 +- 0.022      1.44   138.6     144.8        0.229   100%
cluster <=2ph      1.372 +- 0.022      1.50   136.5     142.7        0.223   100%
cluster <=3ph      1.386 +- 0.018      1.02   135.2     139.7        0.244   100%
GridBased flat     1.346 +- 0.019      1.15   139.2     143.8        0.259   100%

In [5]:

COLOURS = {'eta': '#a0451c', 'cluster 1ph': '#6b6b6b', 'cluster <=2ph': '#7b5aa6',
           'cluster <=3ph': '#4b2a6b', 'GridBased raw': '#4cbcac', 'GridBased flat': '#10695f'}
frequency = np.linspace(0, 400, 2000)
fig, ax = plt.subplots(1, 2, figsize=(12.5, 4.8))
for name, r in results.items():
    ax[0].plot(r['x'], r['z'], lw=1.1, color=COLOURS[name],
               label=f"{name}: $\\sigma$ = {r['sigma']:.3f} $\\pm$ {r['sigma_err']:.3f} $\\mu$m")
    sigma_mm = r['sigma_gauss'] / 1000.0
    ax[1].plot(frequency, np.exp(-2 * np.pi**2 * sigma_mm**2 * frequency**2),
               lw=1.6, color=COLOURS[name], label=f"{name} ({r['mtf50']:.0f} lp/mm)")
ax[0].set_xlim(-9, 9); ax[0].set_ylim(-0.05, 1.35)
ax[0].set_xlabel(r'position across the edge [$\mu$m]')
ax[0].set_ylabel('relative intensity')
ax[0].legend(fontsize=8); ax[0].grid(alpha=.3); ax[0].set_title('oversampled ESF')
ax[1].axvline(200, color='k', ls=':', lw=.8, label='2.5 um Nyquist')
ax[1].axhline(0.5, color='grey', lw=.5)
ax[1].set_xlim(0, 400); ax[1].set_ylim(0, 1.05)
ax[1].set_xlabel('spatial frequency [lp/mm]'); ax[1].set_ylabel('MTF')
ax[1].legend(fontsize=8); ax[1].grid(alpha=.3); ax[1].set_title('MTF from the fitted ESF')
fig.tight_layout()

3. Where the GridBased sub-pixel phase stands¶

In [6]:

before = np.load(GRID / FLAT_RUN / 'subpixel_phase_before.npy')
after = np.load(GRID / FLAT_RUN / 'subpixel_phase_after.npy')

fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))
for ax, histogram, title in ((axes[0], before, 'before flattening'),
                             (axes[1], after, 'after flattening')):
    im = ax.imshow(histogram.T / histogram.mean(), origin='lower', extent=[0, 1, 0, 1],
                   vmin=0.9, vmax=1.1, cmap='viridis')
    ax.set_title(f'flat field, {title}', fontsize=10)
    ax.set_xlabel('sub-pixel x'); ax.set_ylabel('sub-pixel y')
    plt.colorbar(im, ax=ax, fraction=0.046)
for histogram, title in ((before, 'before'), (after, 'after')):
    profile = histogram.sum(axis=1)
    axes[2].plot(np.linspace(0, 1, len(profile)), profile / profile.mean(), lw=1, label=title)
axes[2].axhline(1, color='k', lw=.6); axes[2].set_xlabel('sub-pixel x')
axes[2].set_ylabel('relative density'); axes[2].legend(); axes[2].grid(alpha=.3)
fig.tight_layout()

4. Relative DQE¶

In [7]:

# DQE(f) is proportional to QE * MTF(f)^2.  Every chain sees the same beam, so
# the ratio of detected photons is the ratio of QE and the relative DQE is well
# defined even without the absolute incident flux.
def mtf_at(frequency_lpmm, sigma_gauss_um):
    sigma_mm = sigma_gauss_um / 1000.0
    return float(np.exp(-2 * np.pi**2 * sigma_mm**2 * frequency_lpmm**2))


FREQUENCY = 200.0                      # lp/mm, the Nyquist of the 2.5 um bins
REFERENCE = 'cluster <=3ph'
reference_photons = results[REFERENCE]['photons']
print(f"{'variant':16s} {'QE rel':>8} {'MTF(200)':>9} {'DQE rel':>9}  (relative to {REFERENCE})")
relative_dqe = {}
for name, r in results.items():
    qe = r['photons'] / reference_photons
    mtf200 = r['mtf200_lsf']           # measured from the LSF, as in the poly notebook
    relative_dqe[name] = qe * mtf200**2
    print(f'{name:16s} {qe:8.3f} {mtf200:9.3f} {relative_dqe[name]:9.3f}')
best = max(relative_dqe, key=relative_dqe.get)
print(f"\nhighest relative DQE at {FREQUENCY:.0f} lp/mm: {best} "
      f"({relative_dqe[best]/relative_dqe['eta']:.2f}x eta)")

variant            QE rel  MTF(200)   DQE rel  (relative to cluster <=3ph)
eta                 0.710     0.187     0.025
cluster 1ph         0.710     0.229     0.037
cluster <=2ph       0.929     0.223     0.046
cluster <=3ph       1.000     0.244     0.060
GridBased flat      1.044     0.259     0.070

highest relative DQE at 200 lp/mm: GridBased flat (2.82x eta)

5. Summary¶

In [8]:

summary = {
    'measurement': 'Edge/Flat_2filters_12keV',
    'occupancy_percent': occupancy * 100,                 # flat field, ESF ROI
    'roi_superpixel': list(ROI),
    'bin_um': BIN_UM,
    'bin_width_um': BIN_WIDTH_UM,
    'convention': 'sigma = ODR p2 (Gaussian LSF); MTF50 = 1000 sqrt(ln 2 / 2) / (pi sigma); '
                  'ESF from DeepLearning/src, bin width 0.3 um',
    'grid_model_epoch': grid_summary['epoch'],
    'grid_frames': grid_summary['n_frames'],
    'grid_subpixel_before': grid_summary['subpixel_nonuniformity_before'],
    'grid_subpixel_after': grid_summary['subpixel_nonuniformity_after'],
    'variants': {name: dict(sigma_um=float(r['sigma']), sigma_err_um=float(r['sigma_err']),
                            mtf50_lsf_lpmm=float(r['mtf50_lsf']),
                            mtf200_lsf=float(r['mtf200_lsf']),
                            chi2_ndf=float(r['chi2_ndf']),
                            mtf50_lpmm=float(r['mtf50']), mtf10_lpmm=float(r['mtf10']),
                            rows_kept=float(r['rows_kept']), photons=float(r['photons']))
                 for name, r in results.items()},
    'relative_dqe_200lpmm': relative_dqe,
}
print(json.dumps(summary['variants'], indent=1))
out = GRID / EDGE_RUN / 'knife_edge_comparison.json'
json.dump(summary, open(out, 'w'), indent=1)
print(f'\nwritten to {out}')

{
 "eta": {
  "sigma_um": 1.482743287676068,
  "sigma_err_um": 0.023612628870164847,
  "mtf50_lsf_lpmm": 131.67127170307347,
  "mtf200_lsf": 0.18708166715602986,
  "chi2_ndf": 1.6091412712838244,
  "mtf50_lpmm": 126.38103081415967,
  "mtf10_lpmm": 230.34405458447614,
  "rows_kept": 1.0,
  "photons": 156697901.0
 },
 "cluster 1ph": {
  "sigma_um": 1.3517498054839803,
  "sigma_err_um": 0.02228742348935663,
  "mtf50_lsf_lpmm": 144.76471476919644,
  "mtf200_lsf": 0.22898527725431547,
  "chi2_ndf": 1.4377846589189212,
  "mtf50_lpmm": 138.62818723482948,
  "mtf10_lpmm": 252.66591450992405,
  "rows_kept": 1.0,
  "photons": 156697890.0
 },
 "cluster <=2ph": {
  "sigma_um": 1.3724427091990845,
  "sigma_err_um": 0.021926901151490465,
  "mtf50_lsf_lpmm": 142.69021685073372,
  "mtf200_lsf": 0.2230447974394754,
  "chi2_ndf": 1.5018077194468944,
  "mtf50_lpmm": 136.53803096715998,
  "mtf10_lpmm": 248.85636282080927,
  "rows_kept": 1.0,
  "photons": 205040034.0
 },
 "cluster <=3ph": {
  "sigma_um": 1.3855387111593214,
  "sigma_err_um": 0.017612306263149854,
  "mtf50_lsf_lpmm": 139.7355710728802,
  "mtf200_lsf": 0.2440890142229673,
  "chi2_ndf": 1.0195693691410581,
  "mtf50_lpmm": 135.24748433227265,
  "mtf10_lpmm": 246.50419222530726,
  "rows_kept": 1.0,
  "photons": 220757025.0
 },
 "GridBased flat": {
  "sigma_um": 1.3461043903325456,
  "sigma_err_um": 0.01889890619496615,
  "mtf50_lsf_lpmm": 143.7567164975865,
  "mtf200_lsf": 0.25887598776406284,
  "chi2_ndf": 1.1515016338976902,
  "mtf50_lpmm": 139.2095787481861,
  "mtf10_lpmm": 253.72556782676153,
  "rows_kept": 1.0,
  "photons": 230530601.0
 }
}

written to /mnt/sls_det_storage/moench_data/MLXID/InferenceResults/GridBased/Edge_2filters_12keV/knife_edge_comparison.json