MLXID/DeepLearning
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add funcs for post processing

Browse Source
This commit is contained in:
xie_x1
2026-06-10 11:24:44 +02:00
parent 54e9b4e2df
commit 9a38cdbe8e
1 changed files with 395 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/funcs_postprocess.py
+395
View File
@@ -0,0 +1,395 @@
import numpy as np
from scipy.signal.windows import tukey
from scipy.fft import rfft, rfftfreq
from ROOT import *
from array import array
def generate_presampled_esf(knifeEdge, flatField, roi,
pixel_size_um=2.5,
bin_width_um=0.4,
chi2_threshold=15.0,
residual_tolerance_px=1.0):
"""
基于刀口/刃边法图像提取并构建过采样边缘扩展函数 (Oversampled ESF)。
参数:
- knifeEdge (np.ndarray): 原始刀口图像数据
- flatField (np.ndarray): 平场图像数据 (须与 knifeEdge 尺寸一致)
- roi (tuple): 感兴趣区域 (x_min, x_max, y_min, y_max)
- pixel_size_um (float): 探测器单像素物理尺寸 (微米)
- bin_width_um (float): 输出的均匀重采样网格宽度 (微米，必须满足奈奎斯特采样限制)
- chi2_threshold (float): 初始单行拟合剔除阈值
- residual_tolerance_px (float): 剔除偏离拟合直线异常点的容差 (像素)
返回:
- dict: 包含准备好的 NumPy 数组供 FFT 使用，以及 ROOT TGraphErrors 供直接拟合使用
"""
x_min, x_max, y_min, y_max = roi
# 1. 动态判断边缘方向并定义拟合函数
left_sum = np.average(knifeEdge[y_min:y_max, x_min])
right_sum = np.average(knifeEdge[y_min:y_max, x_max-1])
if left_sum < right_sum:
erf_func = TF1('erf_func', '[0]*0.5*(1 + TMath::Erf((x - [1])/[2])) + [3]')
else:
erf_func = TF1('erf_func', '[0]*0.5*(1 - TMath::Erf((x - [1])/[2])) + [3]')
# 2. 逐行初筛提取边缘亚像素位置
x_midPoints = array('d')
ys = array('d')
skipped_lines = []
for y in range(y_min, y_max):
xs, xserr, zs, zserr = array('d'), array('d'), array('d'), array('d')
for x in range(x_min, x_max):
xs.append(x)
xserr.append(0.0)
n1 = float(knifeEdge[y, x])
n2 = float(flatField[y, x])
safe_n2 = n2 + 1e-6
zs.append(n1 / safe_n2)
zserr.append(np.sqrt(n2**2 * n1 + n1**2 * n2) / (safe_n2**2))
gr = TGraphErrors(len(xs), xs, zs, xserr, zserr)
# 鲁棒初始值估算 (改用排序中位数避免过渡区干扰)
zs_sorted = np.sort(np.array(zs))
p3 = np.median(zs_sorted[:3])
p0 = np.median(zs_sorted[-3:]) - p3
p1 = (x_min + x_max) / 2.0
p2 = 1.8
erf_func.SetParameters(p0, p1, p2, p3)
erf_func.SetParLimits(0, 0.1*p0, 5*p0)
erf_func.SetParLimits(1, x_min, x_max)
erf_func.SetParLimits(2, 0.1, 5.0)
erf_func.SetParLimits(3, -0.5, 5.0)
# Q0W: 静默、不绘制、权重设为1(仅拟合形状)
gr.Fit(erf_func, 'Q0W')
chi2 = erf_func.GetChisquare()
ndf = erf_func.GetNDF()
if ndf == 0 or (chi2 / ndf) > chi2_threshold:
skipped_lines.append(y)
continue
x_midPoints.append(erf_func.GetParameter(1))
ys.append(y)
# 3. 边缘斜线鲁棒拟合与残差剔除
g_xMidpoint_y = TGraph(len(x_midPoints), ys, x_midPoints)
g_xMidpoint_y.Fit('pol1', 'Q0 ROB') # 第一次鲁棒拟合
f_midpoint_rough = g_xMidpoint_y.GetFunction('pol1')
clean_ys, clean_x_midPoints = array('d'), array('d')
for i in range(g_xMidpoint_y.GetN()):
y_val = g_xMidpoint_y.GetPointX(i)
x_val_meas = g_xMidpoint_y.GetPointY(i)
x_val_exp = f_midpoint_rough.Eval(y_val)
if abs(x_val_meas - x_val_exp) < residual_tolerance_px:
clean_ys.append(y_val)
clean_x_midPoints.append(x_val_meas)
# 最终高精度直线拟合
g_clean = TGraph(len(clean_ys), clean_ys, clean_x_midPoints)
g_clean.Fit('pol1', 'Q0 ROB')
f_midpoint = g_clean.GetFunction('pol1')
# 4. 重采样与物理单位转换 (中心强制对齐到 0)
xOversampled, zOversampled, zErrors = [], [], []
for y in range(y_min, y_max):
if y in skipped_lines:
continue
x_mid = f_midpoint.Eval(y)
for x in range(x_min, x_max):
# 将相对中心坐标转换为物理微米 (中心为0)
x_physical_um = (x - x_mid) * pixel_size_um
xOversampled.append(x_physical_um)
k = float(knifeEdge[y, x])
f = float(flatField[y, x])
var_k = max(k, 1.0)
var_f = max(f, 1.0)
f_safe = max(f, 1.0)
zOversampled.append(k / f_safe)
zErrors.append(np.sqrt(f_safe**2 * var_k + k**2 * var_f) / (f_safe**2))
# 按 X 排序
sorted_pairs = sorted(zip(xOversampled, zOversampled, zErrors))
x_arr = np.array([p[0] for p in sorted_pairs])
z_arr = np.array([p[1] for p in sorted_pairs])
z_err_arr = np.array([p[2] for p in sorted_pairs])
# 5. 统一离散分 Bin (用于 FFT 和最终拟合)
x_min_b, x_max_b = x_arr.min(), x_arr.max()
bin_start = np.floor(x_min_b / bin_width_um) * bin_width_um
bin_edges = np.arange(bin_start, x_max_b + bin_width_um * 1.01, bin_width_um)
n_bins = len(bin_edges) - 1
x_binned, z_binned, x_err_binned, z_err_binned = [], [], [], []
for i in range(n_bins):
low, high = bin_edges[i], bin_edges[i+1]
if i == n_bins - 1:
mask = (x_arr >= low) & (x_arr <= high)
else:
mask = (x_arr >= low) & (x_arr < high)
if not np.any(mask): continue
z_vals, z_err_vals, x_vals = z_arr[mask], z_err_arr[mask], x_arr[mask]
z_err_vals = np.clip(z_err_vals, 1e-12, None)
weights = 1.0 / (z_err_vals ** 2)
# 强制网格点在 bin 中心 (FFT 严格要求等距网格)
x_center = low + bin_width_um / 2.0
x_binned.append(x_center)
x_err_binned.append(bin_width_um / np.sqrt(12))
z_binned.append(np.sum(z_vals * weights) / np.sum(weights))
z_err_binned.append(1.0 / np.sqrt(np.sum(weights)))
# 转换为 PyROOT 适用的 array格式
x_root = array('d', x_binned)
z_root = array('d', z_binned)
x_err_root = array('d', x_err_binned)
z_err_root = array('d', z_err_binned)
gr_esf = TGraphErrors(len(x_root), x_root, z_root, x_err_root, z_err_root)
gr_esf.SetTitle('Oversampled ESF; Position [#mum]; Relative Intensity')
gr_esf.SetMarkerStyle(20)
# 返回数据字典
return {
'numpy': {
'x': np.array(x_binned),
'z': np.array(z_binned),
'x_err': np.array(x_err_binned),
'z_err': np.array(z_err_binned)
},
'root_graph': gr_esf,
'edge_slope': f_midpoint.GetParameter(1),
'valid_lines_ratio': len(clean_ys) / (y_max - y_min)
}
def calculate_mtf_iso12233(x_binned, z_binned, pixel_size_um, pad_length=4096):
"""
依据 ISO 12233 标准，从均匀分 bin 的 ESF 计算 MTF，并提取硬件 Nyquist 频率处的 MTF 值。
参数:
x_binned (list/array): 均匀的空间位置网格 (um)
z_binned (list/array): 对应的归一化边缘强度值
pixel_size_um (float): 探测器的物理像素尺寸 (如 2.5)
pad_length (int): FFT 补零长度，必须是2的幂次，用于提高频域分辨率
返回:
freqs_lpmm (array): 空间频率轴 (lp/mm)
mtf (array): 归一化后的 MTF 曲线值
mtf_at_ny (float): 探测器物理 Nyquist 频率处的 MTF 值
nyquist_lpmm (float): 探测器的硬件 Nyquist 频率 (lp/mm)
"""
x_arr = np.array(x_binned)
z_arr = np.array(z_binned)
# 1. 自动获取当前分 bin 的步长 dx
dx_um = x_arr[1] - x_arr[0]
# 2. 数值求导得到 LSF
lsf = np.gradient(z_arr, dx_um)
# 3. LSF 居中
peak_idx = np.argmax(lsf)
half_width = min(peak_idx, len(lsf) - peak_idx - 1)
lsf_centered = lsf[peak_idx - half_width : peak_idx + half_width + 1]
esf_centered = z_arr[peak_idx - half_width : peak_idx + half_width + 1]
# 4. Zero-Padding (补零)
pad_left = (pad_length - len(lsf_centered)) // 2
pad_right = pad_length - len(lsf_centered) - pad_left
lsf_padded = np.pad(lsf_centered, (pad_left, pad_right), 'constant')
# 5. FFT 计算 MTF
mtf_raw = np.abs(rfft(lsf_padded))
mtf = mtf_raw / mtf_raw[0] # 归一化，使得直流分量 MTF(0) = 1.0
# 6. 计算频率轴
freqs_cpum = rfftfreq(pad_length, d=dx_um) # 单位: cycles/um
freqs_lpmm = freqs_cpum * 1000 # 转换为 lp/mm
# 7. 锁定探测器硬件 Nyquist 频率并插值获取对应的 MTF 值
nyquist_lpmm = 1000.0 / (2.0 * pixel_size_um)
mtf_at_ny = float(np.interp(nyquist_lpmm, freqs_lpmm, mtf))
mtf50 = float(np.interp(0.5, mtf[::-1], freqs_lpmm[::-1])) # MTF=0.5 对应的频率
mtf10 = float(np.interp(0.1, mtf[::-1], freqs_lpmm[::-1])) # MTF=0.1 对应的频率
print(f"=== MTF Analysis ===")
print(f"Oversampling Bin Width : {dx_um:.3f} um (Analysis Nyquist: {1000/(2*dx_um):.1f} lp/mm)")
print(f"Detector Super Pixel : {pixel_size_um:.1f} um (Hardware Nyquist: {nyquist_lpmm:.1f} lp/mm)")
print(f"Result -> MTF(Nyquist) : {mtf_at_ny:.2f}")
print(f"Result -> MTF(0.5) : {mtf50:.1f} lp/mm")
print(f"Result -> MTF(0.1) : {mtf10:.1f} lp/mm")
print(f'=====================')
return {
'freqs_lpmm': freqs_lpmm,
'mtf': mtf,
'mtf_at_ny': mtf_at_ny,
'mtf50': mtf50,
'mtf10': mtf10,
'nyquist_lpmm': nyquist_lpmm,
'lsf_centered': lsf_centered,
'esf_centered': esf_centered,
'MTF10_freq_lpmm': mtf10,
'MTF50_freq_lpmm': mtf50,
'lsf_padded': lsf_padded
}
def get_graph_fitting(xPreSampled, zPreSampled, xErrors, zErrors):
gr_oversampled = TGraphErrors(len(xPreSampled), xPreSampled, zPreSampled, xErrors, zErrors)
gr_oversampled.SetMarkerStyle(20)
gr_oversampled.SetTitle(';X Position Aligned [#mum]; Relative Intensity')
gr_oversampled.Draw('APL')
erf_func = TF1('erf_func', '[0]*0.5*(1 + TMath::Erf((x - [1])/[2])) + [3]')
erf_func.SetLineColor(kRed+1)
erf_func.SetParLimits(0, 0.2, 1.2) # amplitude
erf_func.SetParLimits(1, -5*2.5, 5*2.5) # midpoint
erf_func.SetParLimits(2, .6, 3) # sigma
erf_func.SetParLimits(3, -0.2, 0.2) # baseline
gr_oversampled.Fit(erf_func, 'WQ')
fitted_center = erf_func.GetParameter(1)
sigma = erf_func.GetParameter(2)
erf_func.SetRange(fitted_center - 2*sigma, fitted_center + 5*sigma)
gr_oversampled.Fit(erf_func, 'ER')
gr_oversampled.GetYaxis().SetRangeUser(0., 1.4)
gr_oversampled.GetXaxis().SetLimits(fitted_center - 9, fitted_center + 9)
chi2, ndf = erf_func.GetChisquare(), erf_func.GetNDF()
sigma_1ph_um, sigma_1ph_err_um = erf_func.GetParameter(2), erf_func.GetParError(2)
return gr_oversampled, erf_func
from scipy.fft import fft2, fftshift
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
# --- 辅助函数：纯 numpy 实现的高效 Rebin (平局池化) ---
def rebin_2d(arr, bin_factor):
new_shape = (arr.shape[0] // bin_factor, bin_factor,
arr.shape[1] // bin_factor, bin_factor)
return arr.reshape(new_shape).sum(-1).sum(1)
# --- 辅助函数：Logistic 拟合函数 ---
def logistic_func(x, L, k, x0, c):
return L / (1 + np.exp(-k * (x - x0))) + c
def nps_calculator(frames_list, pixel_size_um=2.5, target_cx=120, target_cy=220, target_roi_size=64):
"""
参数:
frames_list: 16张连续时间子图的 list
pixel_size_um: 物理像素尺寸
target_cx, target_cy: ROI中心 (必须确保完全在亮区内)
"""
num_pairs = len(frames_list) // 2
N = target_roi_size
nps_2d_accum = np.zeros((N, N))
diff_rois = [] # 存储所有 ROI 差分图，用于后续计算 NPS(0)
# ---------------------------------------------------------
# 1. 图像相减与 2D FFT
# ---------------------------------------------------------
for i in range(num_pairs):
img1 = frames_list[2*i].astype(np.float64)
img2 = frames_list[2*i + 1].astype(np.float64)
diff_img = img1 - img2
half = N // 2
roi_diff = diff_img[target_cy - half : target_cy + half,
target_cx - half : target_cx + half]
diff_rois.append(roi_diff)
# 去直流分量并 FFT
roi_diff_centered = roi_diff - np.mean(roi_diff)
fft_mag_sq = np.abs(fft2(roi_diff_centered))**2
nps_2d_accum += fft_mag_sq
# 平均 2D 频谱，并使用 "像素归一化" (为了兼容你的 DQE 公式)
# / 2.0 是因为两帧相减方差翻倍
nps_2d_mean = (nps_2d_accum / num_pairs) / (N ** 2) / 2.0
nps_2d_shifted = fftshift(nps_2d_mean)
# ---------------------------------------------------------
# 2. 计算 NPS(0) (Kirsty 的方差外推法)
# ---------------------------------------------------------
max_bin = N // 2
nps0_array = np.zeros(max_bin)
bin_factors = np.arange(1, max_bin + 1)
for b in bin_factors:
dim = N - (N % b)
rebinned_size = dim // b
noise_sum = 0
for roi_diff in diff_rois:
# 提取 4 个象限
q1, q2 = roi_diff[:dim, :dim], roi_diff[-dim:, -dim:]
q3, q4 = roi_diff[:dim, -dim:], roi_diff[-dim:, :dim]
# Rebin 并求方差 (ddof=1)
v1, v2 = np.var(rebin_2d(q1, b), ddof=1), np.var(rebin_2d(q2, b), ddof=1)
v3, v4 = np.var(rebin_2d(q3, b), ddof=1), np.var(rebin_2d(q4, b), ddof=1)
noise_sum += np.mean([v1, v2, v3, v4])
# 平均并除以 b^2 和 2.0
nps0_array[b-1] = (noise_sum / num_pairs) / (b**2) / 2.0
# Logistic 拟合寻找 NPS(0) 渐近线
try:
p0 = [nps0_array[-1] - nps0_array[0], 0.5, max_bin/2, nps0_array[0]]
popt, _ = curve_fit(logistic_func, bin_factors, nps0_array, p0=p0, maxfev=5000)
nps_0_val = popt[0] + popt[3] # L + c (即 x->无穷大 时的极限)
except:
print("Warning: Logistic fit failed. Using average of tail.")
nps_0_val = np.mean(nps0_array[-5:])
# ---------------------------------------------------------
# 3. 径向平均提取 1D NPS (Radial Averaging)
# ---------------------------------------------------------
center = N // 2
Y, X = np.indices((N, N))
R = np.sqrt((X - center)**2 + (Y - center)**2)
R_rounded = np.round(R).astype(int)
nps_1d = np.zeros(center)
counts = np.zeros(center)
for r_val in range(center):
mask = (R_rounded == r_val)
nps_1d[r_val] = np.mean(nps_2d_shifted[mask])
# 替换首位为精准的 NPS(0)
nps_1d[0] = nps_0_val
# ---------------------------------------------------------
# 4. 生成物理频率轴 (lp/mm)
# ---------------------------------------------------------
dx_mm = pixel_size_um / 1000.0
# 由于我们取了 center 长度的半径，频率步长为 1 / (N * dx)
freqs_1d = np.arange(center) / (N * dx_mm)
return freqs_1d, nps_1d, nps_2d_shifted, (bin_factors, nps0_array)
def selfSubPixelNormalization(image):
subPixelDistribution = np.zeros((10, 10))
for i in range(10):
for j in range(10):
subPixelDistribution[i, j] = np.sum(image[i::10, j::10])
subPixelDistribution /= np.mean(subPixelDistribution)
for i in range(10):
for j in range(10):
image[i::10, j::10] /= subPixelDistribution[i, j]
return image, subPixelDistribution