基于opencv和phasecong提取边缘特征点

基于opencv和phasecong提取边缘特征点

文章目录

• 前言
• 一、图像读取后的通道统一
• 二、phasecong函数
• • 2.1 方法介绍
  • 2.2 参数及返回值
• 三、利用opencv的快速检测提取边缘特征点
• 四、复现RIFT_descriptor_no_rotation_invariance函数
• • 4.1 复现RIFT_descriptor_no_rotation_invariance函数
  • 4.2 复现FSC函数
  • 4.3 复现LSM

---

前言

最近在研究关于红外图像和自然光影像匹配的问题，参考到了李加元老师的文章：RIFT: Multi-modal image matching based on radiation-variation insensitive feature transform
.文章源码为matlab，以下整理一下python代码的复现思路

---

一、图像读取后的通道统一 红外图像是单通道，自然管图像为三通道，读取数据后需要将红外影像复制为三通道，使其能够对应到RGB三个通道的特征

if img1.ndim == 2:
    img = np.zeros([len(img1), len(img1[0]), 3], int)
    for i in range(3):
        img[..., i] = img1
    img0 = img1
    img1 = img

if img2.ndim == 2:
    img = np.zeros([len(img2), len(img2[0]), 3], int)
    for i in range(3):
        img[..., i] = img2
    img2 = img

二、phasecong函数 2.1 方法介绍

phasecong是一个对比度不变的边缘和角点检测器，计算图像相位一致性的函数。它最开始是matlab的一个函数（文章中用到的为matlab中的phasecong3），幸有Alistair Muldal大佬实现了python的phasecong,可以直接cv来用。

2.2 参数及返回值

phasecong参数和返回值列表与matlab版本的相同，文章模型只需要获取其返回的M和EO：

参数:
    ------
             
    img             N/A     The input image
    nscale          5       Number of wavelet scales, try values 3-6
    norient         6       Number of filter orientations.
    minWaveLength   3       Wavelength of smallest scale filter.
    mult            2.1     Scaling factor between successive filters.
    sigmaOnf        0.55    Ratio of the standard deviation of the Gaussian
                            describing the log Gabor filter's transfer function
                            in the frequency domain to the filter center
                            frequency.
    k               2.0     No. of standard deviations of the noise energy
                            beyond the mean at which we set the noise threshold
                            point. You may want to vary this up to a value of
                            10 or 20 for noisy images.
    cutOff          0.5     The fractional measure of frequency spread below
                            which phase congruency values get penalized.
    g               10      Controls the 'sharpness' of the transition in the
                            sigmoid function used to weight phase congruency
                            for frequency spread.
    noiseMethod     -1      Parameter specifies method used to determine
                            noise statistics.
                            -1 use median of smallest scale filter responses
                            -2 use mode of smallest scale filter responses
                            >=0 use this value as the fixed noise threshold
返回值:
    ---------
    M       Maximum moment of phase congruency covariance, which can be used as
            a measure of edge strength
    m       Minimum moment of phase congruency covariance, which can be used as
            a measure of corner strength
    ori     Orientation image, in integer degrees (0-180), positive angles
            anti-clockwise.
    ft      Local weighted mean phase angle at every point in the image. A
            value of pi/2 corresponds to a bright line, 0 to a step and -pi/2
            to a dark line.
    PC      A list of phase congruency images (values between 0 and 1), one per
            orientation.
    EO      A list containing the complex-valued convolution results (see
            below)
    T       Calculated noise threshold (can be useful for diagnosing noise
            characteristics of images). once you know this you can then specify
            fixed thresholds and save some computation time.

也可参考matlab PHASECONG3函数解析

m1, __, __, __, __, eo1, __ = phasecong.phasecong(img=img1, nscale=4, norient=6, minWaveLength=3, mult=1.6,
                                                  sigmaOnf=0.75, g=3, k=1)
m2, __, __, __, __, eo2, __ = phasecong.phasecong(img=img2, nscale=4, norient=6, minWaveLength=3, mult=1.6,
                                                  sigmaOnf=0.75, g=3, k=1)

---

三、利用opencv的快速检测提取边缘特征点

使用快速检测算法实现提取点 *** 作
注意：需要先将提取到的特征转为unit8格式，不然会报错

m1, m2 = map(lambda img: (img.astype(np.float) - img.min()) / (img.max() - img.min()), (m1, m2))
cm1 = m1 * 255
cm2 = m2 * 255
fast = cv2.FastFeatureDetector_create(nonmaxSuppression=True, type=cv2.FAST_FEATURE_DETECTOR_TYPE_7_12)
kp1 = fast.detect(np.uint8(cm1), None)
kp2 = fast.detect(np.uint8(cm2), None)

---

四、复现RIFT_descriptor_no_rotation_invariance函数

matlab可以直接 *** 作矩阵，python里面对矩阵的 *** 作主要依赖于numpy的array，array需要从detect()的返回值keypoint对象中提取坐标（但好像python也有matrix，不太了解就没用）

4.1 复现RIFT_descriptor_no_rotation_invariance函数

def RIFT_descriptor_no_rotation_invariance(img, kps, eo, patch_size, s, o):
    KPS = kps.T
    (yim, xim, _) = np.shape(img)
    CS = np.zeros([yim, xim, o], np.float)
    for j in range(o):
        for i in range(s):
            # 将各个scale的变换结果的幅度相加
            CS[..., j] = CS[..., j] + np.abs(np.array(eo[j][i]))
    mim = np.argmax(CS, axis=2)
    des = np.zeros([36 * o, np.size(KPS, 1)])
    kps_to_ignore = np.ones([1, np.size(KPS, 1)], bool)
    for k in range(np.size(KPS, 1)):
        x = round(KPS[0][k])
        y = round(KPS[1][k])
        x1 = max(0, x - math.floor(patch_size / 2))
        y1 = max(0, y - math.floor(patch_size / 2))
        x2 = min(x + math.floor(patch_size / 2), np.size(img, 1))
        y2 = min(y + math.floor(patch_size / 2), np.size(img, 0))

        if y2 - y1 != patch_size or x2 - x1 != patch_size:
            kps_to_ignore[0][i] = 0
            continue

        patch = mim[y1:y2, x1:x2]
        ys, xs = np.size(patch, 0), np.size(patch, 1)
        ns = 6;
        RIFT_des = np.zeros([ns, ns, o])
        for j in range(ns):
            for i in range(ns):
                clip = patch[round((j) * ys / ns):round((j + 1) * ys / ns),
                       round((i) * xs / ns): round((i + 1) * xs / ns)]
                x, __ = np.histogram(clip.T.flatten(), bins=6, range=(0, o), density=False)
                te = RIFT_des[j][i]
                RIFT_des[j][i] = x.reshape(1, 1, len(x))
        RIFT_des = RIFT_des.T.flatten()

        df = np.linalg.norm(RIFT_des)
        if df != 0:
            RIFT_des = RIFT_des / df
        des[:, [k]] = np.expand_dims(RIFT_des, axis=1)
    m = repeat(kps_to_ignore, '1 n -> c n', c=2)
    v = KPS[m]
    KPS_out = v.reshape(2, int(len(v) / 2)).T
    w = repeat(kps_to_ignore, '1 n -> c n', c=len(des))
    z = des[w]
    des_out = z.reshape(len(des), int(len(z) / len(des))).T
    des_out = np.float32(des_out) * 100
    return KPS_out, des_out

4.2 复现FSC函数

（只保留了’affine’部分）

def FSC(cor1, cor2, change_form, error_t):
    (M, N) = np.shape(cor1)
    if (change_form == 'similarity'):
        n = 2
        max_iteration = M * (M - 1) / 2
    elif (change_form == 'affine'):
        n = 3
        max_iteration = M * (M - 1) * (M - 2) / (2 * 3)
    elif (change_form == 'perspective'):
        n = 4
        max_iteration = M * (M - 1) * (M - 2) / (2 * 3)

    if (max_iteration > 10000):
        iterations = 10000
    else:
        iterations = max_iteration

    most_consensus_number = 0
    cor1_new = np.zeros([M, N])
    cor2_new = np.zeros([M, N])

    for i in range(iterations):
        while (True):
            a = np.floor(1 + (M - 1) * np.random.rand(1, n)).astype(np.int)[0]
            cor11 = cor1[a]
            cor22 = cor2[a]
            if n == 2 and (a[0] != a[1]) and sum(cor11[0] != cor11[1]) and sum(cor22[0] != cor22[1]):
                break
            if n == 3 and (a[0] != a[1] and a[0] != a[2] and a[1] != a[2]) and sum(cor11[0] != cor11[1]) and sum(
                    cor11[0] != cor11[2]) and sum(cor11[1] != cor11[2]) and sum(cor22[0] != cor22[1]) and sum(
                cor22[0] != cor22[2]) and sum(cor22[1] != cor22[2]):
                break
            if n == 4 and (
                    a[0] != a[1] and a[0] != a[2] and a[0] != a[3] and a[1] != a[2] and a[1] != a[3] and a[2] != a[
                3]) and sum(cor11[0] != cor11[1]) and sum(cor11[0] != cor11[2]) and sum(cor11[0] != cor11[3]) and sum(
                cor11[1] != cor11[2]) and sum(cor11[1] != cor11[3]) and sum(cor11[2] != cor11[3]) and sum(
                cor22[0] != cor11[1]) and sum(cor22[0] != cor22[2]) and sum(cor22[0] != cor22[3]) and sum(
                cor22[1] != cor22[2]) and sum(cor22[1] != cor22[3]) and sum(cor22[2] != cor22[3]):
                break
        parameters, __ = LSM(cor11, cor22, change_form)
        solution = np.array([[parameters[0], parameters[1], parameters[4]],
                             [parameters[2], parameters[3], parameters[5]],
                             [parameters[6], parameters[7], 1]])
        match1_xy = np.ones([3, len(cor1)])
        match1_xy[:2] = cor1.T

        if change_form == 'affine':
            t_match1_xy = solution.dot(match1_xy)
            match2_xy = np.ones([3, len(cor1)])
            match2_xy[:2] = cor2.T
            diff_match2_xy = t_match1_xy - match2_xy
            diff_match2_xy = np.sqrt(sum(np.power(diff_match2_xy, 2)))
            index_in = np.argwhere(diff_match2_xy < error_t)
            consensus_num = len(index_in)
            index_in = np.squeeze(index_in)

        if consensus_num > most_consensus_number:
            most_consensus_number = consensus_num
            cor1_new = cor1[index_in]
            cor2_new = cor2[index_in]
    unil = cor1_new
    __, IA = np.unique(unil, return_index=True, axis=0)
    IA_new = np.sort(IA)
    cor1_new = cor1_new[IA_new]
    cor2_new = cor2_new[IA_new]
    unil = cor2_new
    __, IA = np.unique(unil, return_index=True, axis=0)
    IA_new = np.sort(IA)
    cor1_new = cor1_new[IA_new]
    cor2_new = cor2_new[IA_new]

    parameters, rmse = LSM(cor1_new, cor2_new, change_form)
    solution = np.array([[parameters[0], parameters[1], parameters[4]],
                         [parameters[2], parameters[3], parameters[5]],
                         [parameters[6], parameters[7], 1]])
    return solution

4.3 复现LSM

def LSM(match1, match2, change_form):
    A = np.zeros([2 * len(match1), 4])
    for i in range(len(match1)):
        A[2 * i:2 * i + 2] = np.tile(match1[i], (2, 2))
    B = np.array([[1, 1, 0, 0], [0, 0, 1, 1]])
    B = np.tile(B, (len(match1), 1))
    A = A * B
    B = np.array([[1, 0], [0, 1]])
    B = np.tile(B, (len(match1), 1))
    A = np.hstack((A, B))
    b = match2.reshape(1, int(len(match2) * len(match2[0]))).T

    if change_form == "affine":
        Q, R = np.linalg.qr(A)
        parameters = np.zeros([8, 1])
        parameters[:6] = np.linalg.solve(R, np.dot(Q.T, b))
        N = len(match1)
        M = np.array([[parameters[0][0], parameters[1][0]], [parameters[2][0], parameters[3][0]]])
        match1_test_trans = M.dot(match1.T) + np.tile([parameters[4], parameters[5]], (1, N))
        match1_test_trans = match1_test_trans.T
        test = match1_test_trans - match2
        rmse = math.sqrt(sum(sum(np.power(test, 2))) / N)
    return np.squeeze(parameters), rmse