阶段保存-国标 GB/T 17444-2013 盲元检测-python

#用于验证数据及调用国标中盲元检测算法及添加计算三维噪声功能，注：文中算法经过AI的修正。供后续程序调用，及参考

# -*- coding: utf-8 -*- """ 国标 GB/T 17444-2013 盲元检测 - 死像元(d)：响应率 R < 0.5 * R_mean - 过热像元(h)：噪声电压 VN > 2 * VN_mean - 采用国标推荐的迭代剔除收敛算法 - 全向量化实现，维度安全（统一一维扁平索引，最后再映射回二维标记） """ import numpy as np import os import cv2 as cv # ==================== 基础工具函数 ==================== def calculate_pixel_mean(file_path, frame_number, height, width): """计算所有帧对应像素的平均值（一维扁平数组，长度 height*width）""" n_pixel = height * width image_mean = np.zeros(n_pixel, dtype=np.float64) files = [os.path.join(file_path, f) for f in os.listdir(file_path)] for file in files: with open(file, 'rb') as fp: data = np.fromfile(fp, dtype=np.uint16) if data.size != n_pixel: raise ValueError(f"文件 {file} 数据尺寸 {data.size} 与 {n_pixel} 不符") image_mean += data.astype(np.float64) image_mean /= frame_number return image_mean # 一维，长度 height*width def calculate_pixel_V(file_dir_T0, file_dir_T, frame_number, height, width, K_gain): """像元响应电压 V_DS_mean = (V_T - V_T0) / K_gain""" v_t0 = calculate_pixel_mean(file_dir_T0, frame_number, height, width) v_t = calculate_pixel_mean(file_dir_T, frame_number, height, width) return (v_t - v_t0) / float(K_gain) def calculate_P(T, T0, D, AD, L): """辐照功率差值 P（单位一致即可，cm/K 按国标公式）""" sigma = 5.673e-12 ratio = L / D n = 1.0 if ratio > 1 else 0.0 return sigma * (T**4 - T0**4) * AD / (4 * (ratio**2) + n) def calculate_3d_noise(data_all): """三维噪声模型：输入 (帧数, 行高, 列宽)，返回 7 个噪声分量字典""" t, v, h = data_all.shape mu = data_all.mean() mu_tvh = data_all.mean(axis=0) # (v, h) mu_tvh_rows = data_all.mean(axis=1) # (t, h) mu_tvh_cols = data_all.mean(axis=2) # (t, v) def var_mean(x): return np.mean(x**2) noise = {} # 1. 时空噪声（固定像素帧间变化） noise['时空噪声σ_tvh'] = np.sqrt(var_mean(data_all - mu_tvh[None, :, :])) # 2. 随机空间噪声（固定图案噪声） noise['随机空间噪声σ_vh'] = np.sqrt(var_mean(mu_tvh - mu)) # 3. 行时间噪声 noise['行时间噪声σ_tv'] = np.sqrt(var_mean(mu_tvh_rows - mu_tvh_rows.mean(axis=0)[None, :])) # 4. 列时间噪声 noise['列时间噪声σ_th'] = np.sqrt(var_mean(mu_tvh_cols - mu_tvh_cols.mean(axis=0)[None, :])) # 5. 行固定噪声 noise['行固定噪声σ_v'] = np.sqrt(var_mean(mu_tvh_cols - mu_tvh_cols.mean(axis=1, keepdims=True))) # 6. 列固定噪声 noise['列固定噪声σ_h'] = np.sqrt(var_mean(mu_tvh_rows - mu_tvh_rows.mean(axis=1, keepdims=True))) # 7. 帧处理噪声（时空耦合残留） residual = data_all - (mu_tvh_rows[:, None, :] + mu_tvh_cols[:, :, None] - mu) noise['帧处理噪声σ_t'] = np.sqrt(var_mean(residual)) return noise def calculate_noise_data(file_path, frame_number, height, width, K_gain): """计算噪声：返回 (平均噪声电压标量, 逐像元噪声一维数组, 三维噪声字典)""" n_pixel = height * width files = [os.path.join(file_path, f) for f in os.listdir(file_path)] data_all = np.zeros((frame_number, height, width), dtype=np.float64) for idx, file in enumerate(files): with open(file, 'rb') as fp: data = np.fromfile(fp, dtype=np.uint16) if data.size != n_pixel: raise ValueError(f"文件 {file} 数据尺寸 {data.size} 与 {n_pixel} 不符") data_all[idx] = data.reshape(height, width).astype(np.float64) # 逐像元帧间标准差（样本标准差 ddof=1） noise_img = data_all.std(axis=0, ddof=1) # (height, width) noise_mean = noise_img.sum() / n_pixel # 平均噪声电压 # 限位：噪声过低时钳位到 0.2，保证数值稳定 noise_img_safe = np.where(noise_img <= 0.2, 0.2, noise_img) noise_result = calculate_3d_noise(data_all) return float(noise_mean), noise_img_safe.ravel(), noise_result # 一维返回 # ==================== 国标迭代盲元检测主算法 ==================== def calculate_recommendation_algorithm( file_dir_T0, file_dir_T, frame_number, height, width, K_gain, T, T0, D, AD, L, save_mask_path=None): """ 按国标推荐方法计算死像元(d)与过热像元(h) 返回: d_list(一维索引数组), h_list(一维索引数组), flag(二维掩膜 height x width) """ n_pixel = height * width # 1. 像元响应电压 & 响应率 V_DS_mean = calculate_pixel_V(file_dir_T0, file_dir_T, frame_number, height, width, K_gain) P = calculate_P(T, T0, D, AD, L) R = (V_DS_mean / P).ravel() if P != 0 else np.zeros(n_pixel) # 2. 噪声数据（一维） VN_mean, VN, noise3D = calculate_noise_data(file_dir_T0, frame_number, height, width, K_gain) VN = VN.ravel() R = R.ravel() # 有效像元掩膜：True=仍参与统计的有效像元 valid = np.ones(n_pixel, dtype=bool) # 初始统计（基于全部像元） R_mean = R.mean() Vn_mean = VN.mean() d_set = set() h_set = set() converged = False for _ in range(100): # 迭代上限，防止异常不收敛 # 当前有效索引 cur = np.where(valid)[0] if cur.size == 0: break R_mean = R[cur].mean() Vn_mean = VN[cur].mean() det_d = 0 for idx in cur: if R[idx] < 0.5 * R_mean: d_set.add(idx) valid[idx] = False det_d += 1 det_h = 0 for idx in np.where(valid)[0]: if VN[idx] > 2.0 * Vn_mean: h_set.add(idx) valid[idx] = False det_h += 1 d_cnt = len(d_set) h_cnt = len(h_set) # 国标收敛判据：本轮新增占比均 <= 0.1% ratio_d = (det_d / d_cnt) if d_cnt > 0 else 0.0 ratio_h = (det_h / h_cnt) if h_cnt > 0 else 0.0 if det_d == 0 and det_h == 0: converged = True break if ratio_d <= 0.001 and ratio_h <= 0.001: converged = True break d_list = np.array(sorted(d_set), dtype=np.int64) h_list = np.array(sorted(h_set), dtype=np.int64) # 最终基于收敛后的有效集重算一次均值（供返回/日志） cur = np.where(valid)[0] R_mean_final = float(R[cur].mean()) if cur.size else 0.0 Vn_mean_final = float(VN[cur].mean()) if cur.size else 0.0 # 盲元标记掩膜（二维） flag = np.zeros((height, width), dtype=np.uint8) for idx in d_list: i = idx // width j = idx % width flag[i, j] = 1 for idx in h_list: i = idx // width j = idx % width flag[i, j] = 1 if save_mask_path is not None: cv.imwrite(save_mask_path, flag) print(f"[INFO] 迭代{'收敛' if converged else '达上限'}") print(f"[INFO] 死像元 d = {len(d_list)}, 过热像元 h = {len(h_list)}, " f"盲元总数 = {len(d_list) + len(h_list)}") print(f"[INFO] 最终 R_mean = {R_mean_final:.4f}, VN_mean = {Vn_mean_final:.4f}") return d_list, h_list, flag,noise3D # ==================== 主入口（参数按需修改）==================== if __name__ == "__main__": width = 640 height = 512 frame_number = 100 K_gain = 1.0 file_dir_T = r'D:\CorrectionPic\T' file_dir_T0 = r'D:\CorrectionPic\T0' T = 308.0 T0 = 293.0 D = 4.0 AD = 5.0 L = 20.0 save_mask = r'D:\CorrectionPic\deadPixMask.png' d_list, h_list, flag, noise3D = calculate_recommendation_algorithm( file_dir_T0, file_dir_T, frame_number, height, width, K_gain, T, T0, D, AD, L, save_mask_path=save_mask) print(f"[INFO] 死像元坐标数: {len(d_list)}, 过热像元坐标数: {len(h_list)}") if len(d_list) > 0: print(f"[INFO] 死像元示例(前10): {d_list[:10].tolist()}") if len(h_list) > 0: print(f"[INFO] 过热像元示例(前10): {h_list[:10].tolist()}") print(noise3D)

运行结果：

[INFO] 迭代收敛 [INFO] 死像元 d = 628, 过热像元 h = 31, 盲元总数 = 659 [INFO] 最终 R_mean = 10364630.6662, VN_mean = 9.5830 [INFO] 死像元坐标数: 628, 过热像元坐标数: 31 [INFO] 死像元示例(前10): [611, 629, 633, 636, 637, 638, 639, 640, 641, 1273] [INFO] 过热像元示例(前10): [0, 1, 2, 3, 628, 11815, 82699, 83072, 85384, 93913] {'时空噪声σ_tvh': 9.78523670323088, '随机空间噪声σ_vh': 191.8465215825749, '行时间噪声σ_tv': 1.2903772476809168, '列时间噪声σ_th': 2.233019690143643, '行固定噪声σ_v': 78.1860079985119, '列固定噪声σ_h': 116.61049681874424, '帧处理噪声σ_t': 131.10979706190486}

数据源为BIN文件格式的图像数据，帧数为100帧