阶段保存-国标 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, dtypenp.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, dtypenp.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(axis0) # (v, h) mu_tvh_rows data_all.mean(axis1) # (t, h) mu_tvh_cols data_all.mean(axis2) # (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(axis0)[None, :])) # 4. 列时间噪声 noise[列时间噪声σ_th] np.sqrt(var_mean(mu_tvh_cols - mu_tvh_cols.mean(axis0)[None, :])) # 5. 行固定噪声 noise[行固定噪声σ_v] np.sqrt(var_mean(mu_tvh_cols - mu_tvh_cols.mean(axis1, keepdimsTrue))) # 6. 列固定噪声 noise[列固定噪声σ_h] np.sqrt(var_mean(mu_tvh_rows - mu_tvh_rows.mean(axis1, keepdimsTrue))) # 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), dtypenp.float64) for idx, file in enumerate(files): with open(file, rb) as fp: data np.fromfile(fp, dtypenp.uint16) if data.size ! n_pixel: raise ValueError(f文件 {file} 数据尺寸 {data.size} 与 {n_pixel} 不符) data_all[idx] data.reshape(height, width).astype(np.float64) # 逐像元帧间标准差样本标准差 ddof1 noise_img data_all.std(axis0, ddof1) # (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_pathNone): 按国标推荐方法计算死像元(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, dtypebool) # 初始统计基于全部像元 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), dtypenp.int64) h_list np.array(sorted(h_set), dtypenp.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), dtypenp.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 rD:\CorrectionPic\T file_dir_T0 rD:\CorrectionPic\T0 T 308.0 T0 293.0 D 4.0 AD 5.0 L 20.0 save_mask rD:\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_pathsave_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帧