红外小目标检测算法QDCT复现（python）

引言

最近在专研红外弱小目标检测算法，QDCT发表于2019年，算法也比较经典。由于作者只发布了编译后的MATLAB代码，于是决定动手复现一下。话不多说，直接上干货。

1. 算法框架

2. 特征选择

先从4个独立分量的特征开始。涉及到算法总共有3个，分别是： steerable filter, Kurtosis, Motion

2.1 steerable filter

若多元函数 $f$在点 $P_0(x_1^0,x_2^0,...,x_n^0,)$ 存在对所有自变量的偏导数，则称向量$gradf$为函数 $f$在点 的梯度，记向量的长度（或模）为
$$gradf=(f_{x_1}(P_0),f_{x_2}(P_0),...,f_{x_n}(P_0))$$
$$|gradf|=\sqrt{[f_{x_1}(P_0){]}^2+[f_{x_2}(P_0){]}^2+...+[f_{x_n}(P_0){]}^2}$$
对于函数来说，变化最剧烈的方向即函数在该点的梯度方向，对于一幅灰度图$Im(x,y)$来说，我们可以将其看成是一个二元函数（$p$为像素值，$Im$为原图）
$$p=f(rows,cols)=Im(rows,cols)$$
边缘是像素值变化相对剧烈的点，且该点的梯度方向是像素值变化最剧烈的方向，我们将与该点梯度方向垂直的方向称为该点的边缘方向。求取图像边缘的过程其实就可以等价于求导，然后筛选导数绝对值较大的位置。图像作为一个二元函数，求导准确来说是求方向导数，沿着不同方向求导，得到的值也不一样。根据多元函数的知识我们可以知道，方向导数其实就是该点梯度与该方向方向余弦的点积（投影）。
$$v=gradf\ast \alpha$$
所以当求导方向与该点梯度平行时，最能突出该点的变化，如果垂直，则完全忽视了该点的变化。这样的性质在边缘提取中可以用来特意的突出某一方向的边缘，是设计steerable filter 的基础.

假设输入图像为$Im(x,y)$,关注的方向角度为$\theta$，则$\alpha =PI/2+\theta$, steerable filter的方向余弦为：$(cos(\alpha ),sin(\alpha ))$.
$$\frac{dIm}{dx}=conv2(Im,g_x)$$
$$\frac{dIm}{dy}=conv2(Im,g_y)$$
其中$g_x，g_y$都是沿$x$方向和$y$方向的卷积核，上式给出的是利用卷积求图像偏导的公式，$g_x，g_y$可以是一阶的高斯导数离散化之后的模板，也可以是诸如sobel沿x和y方向的两个卷积核。因此steerable filter的输出为：
$$g_{\alpha }=cos(\alpha )\ast g_x+sin(\alpha )\ast g_y$$
文章中，作者选了$0^o+90^o$和$45^o+135^o$两个特征。
参考代码如下：

def creat_gauss_kernel(kernel_size=3, sigma=1, k=1):
    if sigma == 0:
        sigma = ((kernel_size - 1) * 0.5 - 1) * 0.3 + 0.8
    X = np.linspace(-k, k, kernel_size)
    Y = np.linspace(-k, k, kernel_size)
    x, y = np.meshgrid(X, Y)
    x0 = 0
    y0 = 0
    gauss = 1 / (2 * np.pi * sigma ** 2) * np.exp(- ((x - x0) ** 2 + (y - y0) ** 2) / (2 * sigma ** 2))
    return gauss
 
def get_steerablefilter(Im_in, angle, sigma,debug= 0):
    # 利用高斯一节差分，求图像的可调滤波器结果
    # 输入参数：
    # Im_in : input image, gray
    # angle : unit/rad  range: 0~ 2*pi
    # output paramenters:
    # Im_out image after steerablefilte
    if(Im_in.dtype== 'uint8'):
        Im_in = Im_in.astype('float')
    gauss = creat_gauss_kernel(3, sigma, 1)
    nn = np.zeros_like(gauss)
    nn[:, 0] = [1 / sigma **2, 1 / sigma ** 2, 1 / sigma**2]
    nn[:, 2] = [-1 / sigma **2, -1 / sigma **2, -1 / sigma ** 2]
    gx = np.multiply(nn, gauss)  # 矩阵点乘
    gy = np.multiply(nn.T, gauss)
    # 任意线性滤波器  把图像和卷积核进行卷积  参数2为卷积深度  -1和原图像一样
    dst_x = cv2.filter2D(Im_in, -1, gx)
    dst_y = cv2.filter2D(Im_in, -1, gy)
    Im_out = dst_x * np.cos(angle) + dst_y * np.sin(angle)
    if(debug):
        if(np.max(dst_x)-np.min(dst_x))>1e-3:
            dst_x_save = (dst_x- np.min(dst_x))/(np.max(dst_x)-np.min(dst_x))*255
            dst_x_save = dst_x_save.astype('uint8')
        else:
            dst_x_save = np.zeros_like(Im_in,dtype='uint8')

        if (np.max(dst_y) - np.min(dst_y)) > 1e-3:
            dst_y_save = (dst_y - np.min(dst_y)) / (np.max(dst_y) - np.min(dst_y)) * 255
            dst_y_save = dst_y_save.astype('uint8')
        else:
            dst_y_save = np.zeros_like(Im_in, dtype='uint8')

        if (np.max(Im_out) - np.min(Im_out)) > 1e-3:
            Im_out_save = (Im_out - np.min(Im_out)) / (np.max(Im_out) - np.min(Im_out)) * 255
            Im_out_save = Im_out_save.astype('uint8')
        else:
            Im_out_save = np.zeros_like(Im_in, dtype='uint8')

        cv2.imwrite('sf_dx'+str(int(angle/math.pi*180))+'.png',dst_x_save,[cv2.IMWRITE_PNG_COMPRESSION, 0])
        cv2.imwrite('sf_dy'+str(int(angle/math.pi*180))+'.png', dst_y_save, [cv2.IMWRITE_PNG_COMPRESSION, 0])
        cv2.imwrite('sf_out'+str(int(angle/math.pi*180))+'.png', Im_out_save, [cv2.IMWRITE_PNG_COMPRESSION, 0])
    return Im_out

2.2 Kurtosis feature map

Kurtosis 实际上计算的目标与周围背景的差异性特征，计算比较简单。其计算公式如下：
$$\begin{gathered} kur(Im(x,y,t))=K-C=\frac{E(X-\mu {)}^4}{{\delta }^4}-C \\ =\frac{1}{mn}\ast \frac{{\sum }_{i=1}^{mn}[Im(x,y,t)-\mu {]}^4}{{\delta }^4}-C \end{gathered}$$
作者在文中，选取$m=n=5$, $\mu$和$\delta$分别为每个图像区域的均值和方差，$C$作者取3.

参考代码如下：

def get_kurtosis(Im_in, kernel_size=(5, 5)):
    # get the kurtosis feature map of input image Im_in
    pad_width = ((kernel_size[0] // 2, kernel_size[0] // 2), (kernel_size[1] // 2, kernel_size[1] // 2))
    # 计算 Kurtosis Feature Map
    kfm = np.zeros_like(Im_in)
    for i in range(pad_width[0][0], Im_in.shape[0] - pad_width[0][1]):
        for j in range(pad_width[1][0], Im_in.shape[1] - pad_width[1][1]):
            patch = Im_in[i - pad_width[0][0]:i + pad_width[0][1] + 1, j - pad_width[1][0]:j + pad_width[1][1] + 1]
            kfm[i, j] = kurtosis(patch.ravel())
    return kfm

2.3 motion feature map

在运动特征的计算方面，作者采用了motion charge (MC)的方法，方法计算过程如下：
1、计算相邻帧之间的帧间差分
$$Mov(x,y,t)=\{\begin{array}{rlr}0,& & ifIm(x,y,t)=Im(x,y,t-1)\\ 1,& & ifIm(x,y,t)\ne Im(x,y,t-1)\end{array}$$
2、计算运动累计变量：
$$C{h}_{mov}(x,y,t)=\{\begin{array}{rlr}C{h}_{min},& & ifMov(x,y,t)=1\\ Ch,& & ifMov(x,y,t)=0\end{array}$$
其中$Ch=min(C{h}_{mov}(x,y,t-1)+A,C{h}_{max})$，$C{h}_{mov}(x,y,0)=C{h}_{max}=255$, $C{h}_{min}=0$
3、对motion charge的结果取反，以突出运动特征。
$$m(x,y,t)=255-C{h}_{mov}(x,y,t)$$

参考代码如下：

def get_motion_charge(Im_in, Im_in_pre, thr=3, debug =1):
    # get the motion charge feature of image
    # input parameter
    # Im_in ,the image of now
    # Im_in_pre, the pre image
    #  thr : the threshold of two image
    image_d= np.abs(Im_in - Im_in_pre)>thr
    motion = image_d*255
    if debug:
        save_scale_img('motion.png',motion)
    return motion

3. 四象特征向量

四象的特征向量示意图如下：

$$q(x,y,t)=m(x,y,t)+f_1(x,y,t)i+f_2(x,y,t)j+f_3(x,y,t)k$$

def constract_quard_feature(Im_in,Im_in_pre,thr,debug= 1):
    sigma = 1
    f1 = get_steerablefilter(Im_in,0,sigma)+get_steerablefilter(Im_in,math.pi/2,sigma)
    f2 = get_steerablefilter(Im_in, math.pi/4, sigma) + get_steerablefilter(Im_in, math.pi*3 / 4, sigma)
    f3 = get_kurtosis(Im_in)
    f0 = get_motion_charge(Im_in,Im_in_pre,thr)
    feature = np.zeros((Im_in.shape[0],Im_in.shape[1],4),dtype='float32')
    feature[:,:,0] = f0
    feature[:, :, 1] = f1
    feature[:, :, 2] = f2
    feature[:, :, 3] = f3
    if(debug):
        save_scale_img('f1.png',f1)
        save_scale_img('f2.png',f2)
        save_scale_img('f3.png', f3)
        save_scale_img('f0.png', f0)
    return feature

4. QDCT和IQDCT变换

1、将QDCT应用于q (x、y、t），得到频率特征图Q (u、v、t）。给定四元数q（x，y，t）的分辨率为M×N，其中M，N分别表示图像的长度和宽度。q（x，y，t）的QDCT由表示为：
$$\begin{array}{rl}QDCT(q(x,y,t))& =Q(u,v,t)\\ & \\ & ={\alpha }_u^M{\alpha }_v^N\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{\mu }_{Q}q(x,y,t)N(u,v,m,n)\end{array}$$
其中：
$$N（u，v，m，n）=cos[\frac{\pi }{M}(m+\frac{1}{2})\mu ]cos[\frac{\pi }{N}(n+\frac{1}{2})v]$$

$${\mu }_Q=-\sqrt{\frac{1}{3}}i-\sqrt{\frac{1}{3}}j-\sqrt{\frac{1}{3}}k$$

$${\alpha }_{\mu }^M=\{\begin{array}{ll}\sqrt{\frac{2}{M}},& \mu \ne 0\\ \sqrt{\frac{1}{M}},& \mu =0\end{array}$$

$${\alpha }_v^N=\{\begin{array}{ll}\sqrt{\frac{2}{N}},& v\ne 0\\ \sqrt{\frac{1}{N}},& v=0\end{array}$$

def get_multi_quad(a,b):
    out = np.zeros_like(a)
    out[0] = a[0]*b[0]-a[1]*b[1]-a[2]*b[2]-a[3]*b[3]
    out[1] = a[0]*b[1]+a[1]*b[0]
    out[2] = a[0]*b[2]+a[2]*b[0]
    out[3] = a[0]*b[3]+a[3]*b[0]
    return out
def get_qdct(feature,debug = 1):
    M = feature.shape[0]
    N = feature.shape[1]
    q_out = np.zeros_like(feature)
    pi= math.pi
    b = [0, -np.sqrt(1 / 3), -np.sqrt(1 / 3),-np.sqrt(1 / 3)]
    for u in range(M):
        if u== 0:
            alpha_u = np.sqrt(1/M)
        else:
            alpha_u = np.sqrt(2/M)
        for v in range(N):
            if v == 0:
                alpha_v = np.sqrt(1 / N)
            else:
                alpha_v = np.sqrt(2 / N)
            for m in range(M):
                for n in range(N):
                    if debug:
                        print('get_qdct')
                        print([u,v,m,n])
                    a = feature[m,n,:]
                    q_out[u,v,:] = q_out[u,v,:]+alpha_u*alpha_v*np.cos(pi/M*(m+1/2)*u)*np.cos(pi/N*(n+1/2)*v)*get_multi_quad(a,b)
    return q_out

2、使用符号函数提取变换区域的显著区域自四元数的符号函数对应于视觉反应，用于删除杂乱背景。　　　　　　　　　　　　　　　　　　　　
$$Q^{\prime }=sgn(Q)=\{\begin{array}{ll}\frac{x_0}{|Q|}+\frac{x_1}{|Q|}i+\frac{x_2}{|Q|}j+\frac{x_3}{|Q|}k,& |Q|\ne 0\\ 0,& \text{|Q| =0}\end{array}$$
其中，$x0、x1、x2、x3$分别是Q的四个分量。$|Q|是四元数的大小，|Q|=\sqrt{Q\cdot \overline{Q}}。\overline{Q}是Q$的复共轭。

def get_normalize_quad(q_in,debug = 1):
    M = q_in.shape[0]
    N = q_in.shape[1]
    q_mod = np.zeros_like(q_in)
    for m in range(M):
        for n in range(N):
            if debug:
                print('get_normalize_quad')
                print([m,n])
            q_m = np.sqrt(np.sum(np.multiply(q_in[m,n,:],q_in[m,n,:])))
            if(q_m<1e-5):
                q_mod[m,n,:] = [0,0,0,0]
            else:
                q_mod[m, n, :] =q_in[m,n,:]/q_m
    return q_mod

3、应用$IQDCT将Q^{\prime }$转换回空间域。
$$\begin{array}{rl}{q}^{\prime }(x,y,t)& =IQDCT(sgn(QDCT(q(x,y,t))))=IQDCT(Q^{\prime })\\ & \\ & =\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}{\alpha }_u^M{\alpha }_v^N{\mu }_q{Q}^{\prime }(u,v,t)N(u,v,m,n)\end{array}$$

def get_iqdct(feature,debug = 1):
   # 计算输入的四元数的反四元数变换
   M = feature.shape[0]
   N = feature.shape[1]
   q_out = np.zeros_like(feature)
   pi = math.pi
   b = [0, -np.sqrt(1 / 3), -np.sqrt(1 / 3), -np.sqrt(1 / 3)]
   for m in range(M):
       for n in range(N):
           for u in range(M):
               if u == 0:
                   alpha_u = np.sqrt(1 / M)
               else:
                   alpha_u = np.sqrt(2 / M)
               for v in range(N):
                   if v == 0:
                       alpha_v = np.sqrt(1 / N)
                   else:
                       alpha_v = np.sqrt(2 / N)
                   if debug:
                       print('get_idct')
                       print([m,n,u,v])
                   a = feature[u, v, :]
                   q_out[m, n, :] = q_out[m, n, :] + alpha_u * alpha_v * np.cos(pi / M * (m + 1 / 2) * u) * np.cos(
                       pi / N * (n + 1 / 2) * v) * get_multi_quad(a, b)
   return q_out

4、为了得到更平滑的结果，对$IQDCT$的结果进行高斯滤波，使用$S(x,y,t)$得到最终重建的目标检测结果。
$$S(x,y,t)=g(x,y)\ast [q^{\prime }(x,y,t)\odot \overline{q^{\prime }}(x,y,t)]$$
其中$g（x，y）$是一个使用$\sigma =1.5$的高斯平滑滤波器。$\odot$表示元素级乘积。

参考代码如下：

def get_qf_norm(q_in,debug =1):
    #计算输入的四元数的模
    M = q_in.shape[0]
    N = q_in.shape[1]
    q_norm = np.zeros((M,N),dtype='float32')
    for m in range(M):
        for n in range(N):
            if debug:
                print('get_qf_norm')
                print([m,n])
            q_norm[m,n]= np.sqrt(np.sum(np.multiply(q_in[m, n, :], q_in[m, n, :])))
    return q_norm

整个算法的参考代码如下：

def target_detector_qdct(Im_in,Im_in_pre,thr):
    feature = constract_quard_feature(Im_in, Im_in_pre, thr)
    q_in    = get_qdct(feature)
    q_mod  = get_normalize_quad(q_in)
    q_out = get_iqdct(q_mod)
    q_norm = get_qf_norm(q_out)
    gauss = creat_gauss_kernel(5, 1.5, 2)
    s_out = cv2.filter2D(q_norm, -1, gauss)
    save_scale_img('s_out.png',s_out)
    return s_out

5、实验结果对比

为了简化实验，实验过程中只选取了两张连续图片。如下所示：


将f0, f1,f2,f3均归一化后，计算过程的结果如下：（PS：qdct算法的计算量相当大，复现的代码的计算效率很低，因此本次实验为将10241024的图像缩放到256256处理后得到的结果。）

实验效果与原作者给出的效果，差异很大。可能是一些细节的处理与作者的不一致，还需要继续优化。也希望跟大家继续探讨。