快速径向对称变换（Fast Radial Symmetry Transform）梳理及python代码

关于快速径向对称变换（以下简称FRST变换）的基本介绍可以参考该博客，FRST变换主要考察一个像素点的邻域对当前像素点的作用。

原理和过程

原理和过程如下：

1. 定义orientation projection image On（方向投影图On） 和 magnitude projection image （尺寸投影图Mn）并将其初始化为0.
2. 计算输入图像上各点的梯度值，可以通过计算X方向上和Y方向上的梯度值然后叠加来实现。n是在某一像素位置的搜索范围。那么某一点的图像梯度是有指向的（当然也可能有0的情况），那么梯度g指向的点即为积极影响点（positively-affected pixel），背离的方向即为消极影响点（negatively-affected pixel）。
3. 根据以下公式就可以就算积极影响点和消极影响点的位置其中round是四舍五入取整函数。
4. 对每个像素点进行判断，如果是积极点就在方向投影图相应位置加1，反之则减1。尺寸投影图则在相应位置增减梯度的模。
5. 那么半径为n的径向分布就可以定义为下式
   $A_n$就则是二维的高斯滤波算子，其目的则是让特征点更加明显。$F_n$的计算如下
   
   $\alpha$实际上是矩阵$O_n$的幂，$\alpha$越大，对特征点是对称点的要求越高，论文建议我们一般取2。整个变换可以看做为各个尺寸下得到的变换的总和。
   
   在这里文章建议可以采用奇数尺寸，即n=1,3,5…累加，得到的效果与自然数尺寸直接累加的结果相似。与此同时，也可以对梯度值较低，即变换较为平缓的部分略去（可通过设置阈值来实现，文中阈值系数为$\beta$，取0.2-0.4为宜）。这样可以节省一部分的计算量。

python程序

python程序主要来源于github,在此基础上进行了一定的修改，1.以保证即使在不同的搜索半径的情况下，输出的图片尺寸依然大小相等，便于后续的叠加处理；2.卷积计算采用的是文献推荐的8和9.9两个系数，效果更好。另外提醒输入图片一定要处理为灰度图像。

'''
Implementation of fast radial symmetry transform in pure Python using OpenCV and numpy.
Adapted from:
https://github.com/Xonxt/frst
Which is itself adapted from:
Loy, G., & Zelinsky, A. (2002). A fast radial symmetry transform for detecting points of interest. Computer Vision, ECCV 2002.
'''

import cv2
import numpy as np
import matplotlib.pyplot as plt

def gradx(img):
  img = img.astype('int')
  rows, cols = img.shape
  # Use hstack to add back in the columns that were dropped as zeros
  return np.hstack( (np.zeros((rows, 1)), (img[:, 2:] - img[:, :-2])/2.0, np.zeros((rows, 1))) )

def grady(img):
    img = img.astype('int')
    rows, cols = img.shape
    # Use vstack to add back the rows that were dropped as zeros
    return np.vstack( (np.zeros((1, cols)), (img[2:, :] - img[:-2, :])/2.0, np.zeros((1, cols))) )

#Performs fast radial symmetry transform
#img: input image, grayscale
#radii: integer value for radius size in pixels (n in the original paper); also used to size gaussian kernel
#alpha: Strictness of symmetry transform (higher=more strict; 2 is good place to start)
#beta: gradient threshold parameter, float in [0,1]
#stdFactor: Standard deviation factor for gaussian kernel
#mode: BRIGHT, DARK, or BOTH
def frst(img, radii, alpha, beta, stdFactor, mode='BOTH'):
    mode = mode.upper()
    assert mode in ['BRIGHT', 'DARK', 'BOTH']
    dark = (mode == 'DARK' or mode == 'BOTH')
    bright = (mode == 'BRIGHT' or mode == 'BOTH')

    workingDims = tuple((e + 2*radii) for e in img.shape)#1

#Set up output and M and O working matrices
    output = np.zeros(img.shape, np.uint8)
    O_n = np.zeros(workingDims, np.int16)#每次计算开始需要初始化
    M_n = np.zeros(workingDims, np.int16)

#Calculate gradients
    gx = gradx(img)
    gy = grady(img)

  #Find gradient vector magnitude
    gnorms = np.sqrt(np.add(np.multiply(gx, gx), np.multiply(gy, gy)))

  #Use beta to set threshold - speeds up transform significantly
    gthresh = np.amax(gnorms)*beta

  #Find x/y distance to affected pixels
    gpx = np.multiply(np.divide(gx, gnorms, out=np.zeros(gx.shape), where=gnorms != 0), radii).round().astype(int)
    gpy = np.multiply(np.divide(gy, gnorms, out=np.zeros(gy.shape), where=gnorms != 0), radii).round().astype(int)

  #Iterate over all pixels (w/ gradient above threshold)
    for coords, gnorm in np.ndenumerate(gnorms):
        if gnorm > gthresh:
            i, j = coords
      #Positively affected pixel
            if bright:
                ppve = (i+gpx[i,j], j+gpy[i,j])
                O_n[ppve] += 1
                M_n[ppve] += gnorm
      #Negatively affected pixel
            elif dark:
                pnve = (i-gpx[i,j], j-gpy[i,j])
                O_n[pnve] -= 1
                M_n[pnve] -= gnorm
        else:
            continue

    if radii == 1:
        kn = 8
    else:
        kn = 9.9
  #Abs and normalize O matrix
    O_n = np.abs(O_n)
    O_n = O_n / kn

  #Normalize M matrix
    M_max = float(np.amax(np.abs(M_n)))
    M_n = M_n / kn

#Elementwise multiplication
    F_n = np.multiply(np.power(O_n, alpha), M_n)

#Gaussian blur
    kSize = int(np.ceil(radii / 2))
    if kSize % 2 == 0:
        kSize = kSize + 1
    else:
        kSize = kSize
    S = cv2.GaussianBlur(F_n, (kSize, kSize), int(radii * stdFactor))
    return S[radii:-radii, radii:-radii]