diff --git a/src/funcs_postprocess.py b/src/funcs_postprocess.py new file mode 100644 index 0000000..bc443c1 --- /dev/null +++ b/src/funcs_postprocess.py @@ -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 \ No newline at end of file