时域高通滤波算法（THPF）

时域高通滤波算法（THPF）

今天介绍一下图像非均匀校正中的场景校正算法-时域高通滤波算法（Temporal High Pass Filtering, THPF)以及它后续的各种改进版本。空域和频域已经在之前的文章介绍过一些了，时域还没有。图像是二维空间域上的像素，随着时间的延续，每秒25、30、60帧，就成了视频，针对视频进行图像处理，就可以考虑增加一个时间维度，图像帧与帧之间是存在相关性的，结合时域滤波算法可以有效去除图像中的噪声和探测非均匀性问题。

红外焦平面阵列 (Infrared Focal Plane Arrays，IRFPAs) 由于制造工艺和外部因素影响存在响应非均匀 性问题， 即其所构成的成像系统即使在相同辐射输入的条件下，探测单元的输出也不一致，最终在输出图像上呈现出固定图案噪声，从而严重地影响了成像质量和图像信噪比。目前对于非均匀性校正主要分为两类：一类是基于参考辐射源的定标算法，如两点校正法和多点校正法；另一类是基于场景 的校 正 算 法 ，如时域高通滤波法、神经网络法、恒定统计法和卡尔曼滤波法。后者是利用对场景信息的估计来实现非均匀性校正， 不需要暂停成像系统的工作来重复定标，因此是目前的研究热点。

1. 时域高通滤波非均匀性校正算法

时域高通滤波法是一个时间域的统计滤波器，其基本思想是假设非均匀性处在背景杂波之中，且在图像中随着时间变化相对比较缓慢，而人眼感兴趣的目标信号相对于背景杂波在像平面上运动比较快，所以可以通过时域低通滤波的方法将含噪图像在时域上分解成目标信号( 高频分量) 和背景杂波(低频分量) 两部分。将原始信号与低通滤波器的结果取差即可得到高频分量，即目标信号，从而构成一个高通滤波器，其具体公式表述如下：


通过时域不断迭代得到经过校正第n帧的输出图像，实现过程很简单很经典，效果一般，快速移动目标会有“鬼影”残留，所以需要后续不断改进。

2. 空域低通时域高通非均匀性校正算法

由于影响THPF算法收敛速度与其鬼影问题的根本原因在于大量无关的场景信息被参与到非均匀性校正参数的计算过程中，所以只要尽可能将更多的场景信息，特别是强物体从原始图像中排除，再将剩下的部分参与到非均匀性校正参数的计算过程，这样就可以尽可能减小非随机运动与场景中强物体对校正过程的影响，有效减少鬼影效应。根据这个思想，SLPF-NUC 预先采用空域滤波器将输入的图像信号进行分离。由于在空间上，场景中目标信号具有连续相关性，表现为空域低频性，而非均匀性是由各探测器阵列元独立产生，各阵列元对应的输出的非均匀性差异较大，表现为空域高频性，所以应用空域滤波器将原始信号分离成高频( HSF) 和低频( LSF) 两个部分，只利用高频部分参与非均匀性的校正。

最早的算法中采用了空域线性均值滤波器，场景中的边缘信号同时被平均，所以可以引入一个阈值Th，大于阈值时，可认为它是场景的边缘，置零，不参与非均匀性的校正。后续不断升级低通滤波器，SLPF(空域均值滤波高通非均匀性校正算法)->BFTH(空域双边滤波高通非均匀性校正算法)->GFTH(空域引导滤波高通非均匀性校正算法)->NLMTH(空域非局部均值滤波高通非均匀性校正算法)。也可以通过添加运动检测判断、自适应权重系数对算法进行改进，效果越来越好的同时，算法实时硬件实现的难度也在不断提高。

MATLAB代码实现：

clear; 
close all;
clc;

%% 初始化
frameWidth  = 640;
frameHeight = 512;
M = 200;      % 时间参数
image(:,:,:) = zeros(512,640,100); 
image_low(:,:,:) = zeros(512,640,100); 
image_high(:,:,:) = zeros(512,640,100); 
hsize = 5;
H = fspecial('average',hsize);

%% 帧间迭代
for k = 2:100
    %读取序列图像
    fileName = ['.\序列图像\(', num2str(k), ')', '.Raw'];
   image(:,:,k) = double(reshape(uint16(fread(fopen(fileName),'uint16'))
   				,[frameWidth frameHeight])');   
   image_low(:,:,k)=filter2(H,image(:,:,k));
   image_high(:,:,k) = image(:,:,k) - image_low(:,:,k);
   image(:,:,k) = 1/M * image_high(:,:,k) + (1-1/M)*image(:,:,k-1);  
end 

image_x(:,:) = image(:,:,k);

%% 测试图像
filename    = '第100帧.Raw';
fid         = fopen(filename);
image_raw   = fread(fid,'uint16');
image_raw   = uint16(image_raw);
fclose(fid);
data_temp  = reshape(image_raw,[frameWidth frameHeight]);
image_raw   = data_temp';
image_raw = double(image_raw);
image_out = image_raw - image_x;

%% 输出结果图
figure(1);imshow(image_raw,[]);title('原图像');
figure(2);imshow(image_out,[]);title('SLPF效果图');