基于多向梯度局部对比度的红外弱小目标检测算法（MDWCM）解析及代码（2021年IEEE GRSL 论文）

该论文是< IEEE GEOSCIENCE AND REMOTE SENSING LETTERS >期刊2020年录取的一篇较新的一篇论文《Robust Infrared Small Target Detection via Multidirectional Derivative-Based Weighted Contrast Measure》，目前网上下载地址为：https://www.researchgate.net/publication/344529590_Robust_Infrared_Small_Target_Detection_via_Multidirectional_Derivative-Based_Weighted_Contrast_Measure

探测器距离目标较远的时候，由于目标的尺寸小、亮度弱，场景中由存在较大的背景干扰，使得复杂场景下的弱小点目标检测成为计算机视觉领域一项极具挑战的研究任务。虽然弱小目标亮度较弱，但在局部区域仍然是有一定的对比特性，为此，该论文在分析弱小目标空间分布的基础上，设计了一种基于多向梯度特性与局部对比度的红外弱小点目标检测算法，通过Facet模型实现对图像的多梯度方向的分析，构建了多梯度方向置信图，然后设计了一种加权的局部对比度量算法，有效的区分了目标和背景杂波，最后通过阈值分割实现了对目标的有效检测。算法的流程框架如图1所示。

一、多方向梯度特性计算

由于红外弱小目标可看作图像中的孤立奇异点，相比于较为平坦的背景区域（B）、具有特定方向的边缘文理（E）和噪声（PNHB）等干扰，小目标在各个方向上的梯度特性具有自身明显的特点，论文中利用Facet模型来设计图像的多向梯度滤波器，实现对小目标的多方向分析。

Facet小面模型最早是由Robert M.Haralick提出，该模型提供了一种对邻域图像快速而准确的拟合方法。经典的Haralick Facet小面模型认为图像中每个Facet区域R的灰度分布可由一组双三次多项式基函数的线性组合逼近，逼近函数为：

（1）

通常采用最小二乘方法计算加权系数，但是通过公式(1)求得的解是病态的，而且计算量很大。为此，一般考虑采用2个一维离散chebyshev多项式集合构造二维离散正交多项式基以减少运算量，降低计算过程的病态程度：

（2）

将上式表述为，其中。根据基的正交性和上述的最小均方误差，便可以求得加权系数：

（3）

其中。上式表明加权系数可以通过图像与固定模板卷积获得，其中为模板的大小。

为了提高鲁棒性，采用均值方向梯度来表示图像中一点的梯度。一幅图像的Facet小面拟合曲面如公式（2），对于给定方向矢量，那么在点处沿矢量的一阶方向导数为：

（4）

那么x方向和y方向的方向导数可以求得：

（5）

（6）

取，便可以得到红外图像在四个方向的梯度分布图。

二、MDWCM构建

对于高斯分布的弱小目标，如图2所示。可以看出，在四个方向子带上，弱小目标的有着极具明显的分布特征。

为此，该论文利用了弱小目标的空间分布特性，对弱小目标的局部邻域进行了空间划分，如图3所示。

可以看出，局部的空间划分可以更好的捕获弱小目标的梯度特性。

在空间划分的基础上，作者构建了一种基于梯度特性的局部对比度来增强目标：

为了更加有效的增强目标区域，抑制背景区域。作者并没有使用常用的熵来进行加权，而是利用了一种更加激进的标准差来加权：

最后，就得到了作者构建的MDWCM算子的模型：

对四个梯度方向自带进行乘积操作，从而来进一步更强目标信息：

三、目标分割与检测

按照常规的自适应预测分割算法对最后的显著图进行分割，从而来判断目标。

四、仿真试验

作者在12个图像序列种进行了测试。并对比了一些经典的算法（LCM,MFMM,MPCM,WLDM,DECM,IPT），表面了算法了在某型情况下，具有较好的性能。

算法的计算量较小，实时性也较为满意。

五、总结

该论文是是2020年录用的很新的一篇论文，作者利用了梯度特性+局部对比特性，设计了一种新颖的对比度算子，并具有相对满意的结果。但是基于梯度类的弱小目标检测，在一些复杂的场景下，特别是地面场景和高亮区域的干扰，结果并不能让人满意。

附录：自己编写的代码

%% 输入灰度图像，输出三维图像
function Ct = main_ours(f)

f= double(f);
[M,N] = size(f);

%% 多梯度方向特征图
W2 = [-2,-2,-2,-2,-2;
      -1,-1,-1,-1,-1;
      0,0,0,0,0;
      1,1,1,1,1;
      2,3,4,5,6]./50;
W7 = [-1,-1,-1,-1,-1;
      2,2,2,2,2;
      0,0,0,0,0;
      -2,-2,-2,-2,-2;
      1,1,1,1,1]./60;
W8 = [-4,-2,0,2,4;
    2,1,0,-1,-2;
    4,2,0,-2,-4;
    -2,-1,0,1,2;
    4,2,0,-2,-4]./140;
W3 = W2';
W9 = W8';
W10 = W7';

dfr = imfilter(f,W2,'replicate')-17*imfilter(f,W7,'replicate')/5-2*imfilter(f,W9,'replicate');
dfc = imfilter(f,W3,'replicate')-17*imfilter(f,W10,'replicate')/5-2*imfilter(f,W8,'replicate');

ind = 0;
dfg = zeros(M,N,4);
for t=0:45:359
    ind = ind+1;
    dfg(:,:,ind) = dfr*sind(t)+dfc*cosd(t);
end

%% 0度的滤波器设置
H11 = [0,1,1,1,1;
      0,1,1,1,1;
      0,1,1,1,1];
H12 = [0,0,0,0,0,1,0,0,0,0,-1];  
f11 = imfilter(dfg(:,:,1),H11./12);%均值滤波器
f12 = imfilter(dfg(:,:,1),H12);%差值滤波器
f13 = stdfilt(f12,H11);%差值的标准差
%% 90度的滤波器设置，用于快速运算
H31 = [0,1,1,1,1;
      0,1,1,1,1;
      0,1,1,1,1]';
H32 = [0,0,0,0,0,1,0,0,0,0,-1]';  

f31 = imfilter(dfg(:,:,3),H31./12);%均值滤波器
f32 = imfilter(dfg(:,:,3),H32);%差值滤波器
f33 = stdfilt(f32,H31);%差值的标准差
%% 45度的滤波器设置，用于快速运算
H21 = [1,1,1;
      1,1,1;
      1,1,1];
H22 = [1,1,1,1,1;
      1,1,1,1,1;
      1,1,0,0,0;
      1,1,0,0,0;
      1,1,0,0,0];
H23 = [0,0,0,1,1;
      0,0,0,1,1;
      0,0,0,1,1;
      1,1,1,1,1;
      1,1,1,1,1];
H24 = [1,1,1,1,1;
      1,1,1,1,1;
      0,0,0,0,0];  
H25 = [1,1,0;
      1,1,0;
      1,1,0]; 
temp21 = zeros(M,N);
temp21(1:end-7,1:end-4) = dfg(8:end,5:end,2);
temp22 = zeros(M,N);
temp22(1:end-2,1:end-7) = dfg(3:end,8:end,2);
temp23 = zeros(M,N);
temp23(1:end-2,1:end-2) = dfg(3:end,3:end,2);
f21 = imfilter(dfg(:,:,2),H21./9);%均值滤波器
f22 = imfilter(dfg(:,:,2),H22./16);%均值滤波器
f23 = imfilter(dfg(:,:,2),H23./16);%均值滤波器
f24 = stdfilt(dfg(:,:,2)-temp21,H24);%差值的标准差
f25 = stdfilt(dfg(:,:,2)-temp22,H25);%差值的标准差
f26 = stdfilt(dfg(:,:,2)-temp23,H21);%差值的标准差
%% 135度的滤波器设置，用于快速运算
H41 = [1,1,1;
      1,1,1;
      1,1,1];
H42 = [1,1,1,1,1;
      1,1,1,1,1;
      0,0,0,1,1;
      0,0,0,1,1;
      0,0,0,1,1];
H43 = [1,1,0,0,0;
      1,1,0,0,0;
      1,1,0,0,0;
      1,1,1,1,1;
      1,1,1,1,1];
H44 = [1,1,1,1,1;
      1,1,1,1,1;
      0,0,0,0,0];  
H45 = [0,1,1;
      0,1,1;
      0,1,1]; 
temp41 = zeros(M,N);
temp41(1:end-7,5:end) = dfg(8:end,1:end-4,4);
temp42 = zeros(M,N);
temp42(1:end-2,8:end) = dfg(3:end,1:end-7,4);
temp43 = zeros(M,N);
temp43(1:end-2,3:end) = dfg(3:end,1:end-2,4);
f41 = imfilter(dfg(:,:,4),H41./9);%均值滤波器
f42 = imfilter(dfg(:,:,4),H42./16);%均值滤波器
f43 = imfilter(dfg(:,:,4),H43./16);%均值滤波器
f44 = stdfilt(dfg(:,:,4)-temp41,H44);%差值的标准差
f45 = stdfilt(dfg(:,:,4)-temp42,H45);%差值的标准差
f46 = stdfilt(dfg(:,:,4)-temp43,H41);%差值的标准差

%% 多方向检测结果
s = 2;
C1 = zeros(M,N);
C2 = zeros(M,N);
C3 = zeros(M,N);
C4 = zeros(M,N);

for i=2*s+1+3:M-2*s-3
    for j=2*s+1+3:N-2*s-3 
        dis1 = min(f11(i-1,j-3)-f11(i-1,j+2),f11(i+1,j-3)-f11(i+1,j+2));
        sig1 =  f13(i-1,j-3)*f13(i+1,j-3);
        C1(i,j) = sig1*dis1;
        
        dis3 = min(f31(i-3,j-1)-f31(i+2,j-1),f31(i-3,j+1)-f31(i+2,j+1));
        sig3 =  f33(i-3,j-1)*f33(i-3,j+1);
        C3(i,j) = sig3*dis3;
        
        dis2 = min(f22(i-2,j-2)-f21(i+1,j+1),f21(i-1,j-1)-f23(i+2,j+2));
        sig2 = (max(f24(i-3,j),f25(i-1,j-3)))*f26(i-1,j-1); 
        C2(i,j) = sig2*dis2;
       
        dis4 = min(f42(i-2,j+2)-f41(i+1,j-1),f41(i-1,j+1)-f43(i+2,j-2));
        sig4 = (max(f44(i-3,j+2),f45(i-1,j+3)))*f46(i-1,j+1);      
        C4(i,j) = sig4*dis4;
    end
end

C1 = max(C1,0);
C1 = C1./max(C1(:));
C3 = max(C3,0);
C3 = C3./max(C3(:));
C2 = max(C2,0);
C2 = C2./max(C2(:));
C4 = max(C4,0);
C4 = C4./max(C4(:));

Ct = zeros(M,N);


Ct=C1.*C2.*C3.*C4;

figure(101),
mesh(Ct)
title('Ours-3D')

end

版权声明：本文为CSDN博主「无聊计算机视觉程序员」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_34148844/article/details/116780764