一、前言

RIFT: Multi-modal Image Matching Based on Radiation-invariant Feature Transform （2020）

在图像配准和特征匹配的任务中，传统的算法如 SIFT（尺度不变特征变换） 主要依赖于图像的强度或梯度信息来检测和描述特征点。然而，这些方法在处理非线性辐射畸变（Nonlinear Radiation Distortion, NRD） 的多模态图像时表现不佳。为此，研究者在2020年提出了一种全新的特征匹配算法——RIFT（Radiation-Invariant Feature Transform），以应对在多模态图像中出现的辐射失配问题。

RIFT方法具备以下几个核心创新点：

• 相位一致性特征检测：与传统强度或梯度为基础的特征检测方法不同，RIFT使用相位一致性（Phase Congruency, PC）来检测特征点。PC是一种对辐射变化不敏感的特征，因此在多模态图像下具有更高的稳定性。RIFT结合角点与边缘信息，提高了特征点的数量和重复性。
• 最大索引图（MIM）描述子：RIFT引入最大索引图（Maximum Index Map,MIM），利用Log-Gabor滤波器在不同尺度和方向上的响应来生成描述子，增强了对NRD的鲁棒性。
• 旋转不变性分析：RIFT对MIM值在旋转情况下的变化进行了分析，确保了算法在图像旋转情况下的匹配能力。

通过在六类多模态图像数据集（如光学-红外、SAR-光学、深度图等）上的实验证明，RIFT相比SIFT、SAR-SIFT在匹配准确性与鲁棒性方面表现更优，尤其适用于存在非线性辐射差异的图像匹配任务。

二、代码实现

本文以下部分提供了RIFT算法在Python中的实现，涵盖特征提取、关键点检测、描述子构建、特征匹配及几何变换估计等完整流程。Github仓库：https://github.com/baymin1/RIFT-demo

2.1 整体流程

下面是RIFT算法的完整流程，包括：

• 使用 phasecong() 提取相位一致性特征
• 使用 FAST 角点检测器获取关键点
• 利用RIFT_descriptor_no_rotation_invariance() 构建特征描述子
• 使用FlannBasedMatcher 进行匹配
• 实现自定义RANSAC (fsc) + LSM 最小二乘拟合估计几何变换
• 针对输入图像进行角度迭代，选取匹配误差最小的旋转角作为最优姿态

2.2 核心代码解读

RIFT函数包含特征提取到变换估计的完整流程：

import cv2
import numpy as np
from FSC import fsc
from phasepack import phasecong
from RIFT_Descriptor_No_Rotation_Invariance import RIFT_descriptor_no_rotation_invariance


def rift(img1, img2):
    # 计算图像相位一致性：m1：相位一致协方差的最大矩
    # eo1 和 eo2 包含了通过复值卷积提取出的图像局部相位信息，这些信息在计算图像的特征描述子时被用来增强描述子的鲁棒性和表达能力。
    m1, __, __, __, __, eo1, __ = phasecong(img=img1, nscale=4, norient=6, minWaveLength=3, mult=1.6,
                                            sigmaOnf=0.75, g=3, k=1)
    m2, __, __, __, __, eo2, __ = phasecong(img=img2, nscale=4, norient=6, minWaveLength=3, mult=1.6,
                                            sigmaOnf=0.75, g=3, k=1)

    # 归一化
    m1, m2 = map(lambda img: (img.astype(np.float64) - 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)  # keypoint  fast检测参数1应该为一个单通道图像，故转为unit8格式
    kp2 = fast.detect(np.uint8(cm2), None)

    # 对关键点根据其响应值由大到小进行排序
    kp1 = sorted(kp1, key=lambda kp: kp.response, reverse=True)
    kp2 = sorted(kp2, key=lambda kp: kp.response, reverse=True)

    # 选择前5000个关键点
    top_kp1 = kp1[:min(5000, len(kp1))]
    top_kp2 = kp2[:min(5000, len(kp2))]

    # # 绘制关键点
    # img1_keypoints = cv2.drawKeypoints(img1, top_kp1, None, color=(0, 255, 0))
    # img2_keypoints = cv2.drawKeypoints(img2, top_kp2, None, color=(0, 255, 0))
    # cv2.imshow('Keypoints', img1_keypoints)
    # cv2.imshow('Keypoints2', img2_keypoints)
    # cv2.waitKey(0)
    # cv2.destroyAllWindows()

    # 得到关键点的坐标矩阵
    m1_point = np.array([kp.pt for kp in top_kp1])
    m2_point = np.array([kp.pt for kp in top_kp2])

    # RIFT描述子对特征点进行描述
    kps1, des1 = RIFT_descriptor_no_rotation_invariance(img1, m1_point, eo1, 96, 4, 6)
    kps2, des2 = RIFT_descriptor_no_rotation_invariance(img2, m2_point, eo2, 96, 4, 6)

    # 生成匹配，获取匹配点坐标序列
    bf = cv2.FlannBasedMatcher()
    # matches里面有两幅图的特征点索引和相似度
    matches = bf.match(des1, des2)
    match_point1 = np.zeros([len(matches), 2], int)
    match_point2 = np.zeros([len(matches), 2], int)

    # 根据关键点的索引找到其在kps1中的坐标
    for m in range(len(matches)):
        match_point1[m] = kps1[matches[m].queryIdx]
        match_point2[m] = kps2[matches[m].trainIdx]

    # 去重，得到去重后的坐标数组，IA是在matches中点对的索引。
    match_point2, IA = np.unique(match_point2, return_index=True, axis=0)
    match_point1 = match_point1[IA]

    #求变换矩阵
    transform = fsc(match_point1, match_point2, 'affine', 2)

    # 计算误差
    Y_ = np.ones([3, len(match_point1)])
    Y_[:2] = match_point1.T
    Y_ = transform.dot(Y_)  # 得到图1内点进行变换后的坐标矩阵

    Y_[0] = Y_[0] / Y_[2]
    Y_[1] = Y_[1] / Y_[2]

    threshold = 50
    error = np.sqrt(sum(np.power((Y_[0:2] - match_point2.T), 2)))  # 每个点对的误差
    inliersIndex = np.squeeze(np.argwhere(error < threshold))

    # 将所有小于门限值的点对误差做均值。
    filtered_errors = error[error < threshold]
    mean_error = np.mean(filtered_errors)

    return mean_error, inliersIndex, match_point1, match_point2

关键处理包括：

• 相位一致性（PC）图像归一化处理；
• FAST特征点提取并排序，筛选Top-N响应关键点；
• 构建描述子并匹配；
• 使用FSC实现基于一致性验证的变换参数估计；
• 计算匹配误差，提取内点。

FSC函数是对经典RANSAC的改进，结合最小二乘拟合模型，对关键点对进行一致性检验，输出仿射或透视变换矩阵

import numpy as np
from LSM import lsm


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 = int(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 (cor11[0] != cor11[1]) and (cor22[0] != cor22[1]):
                break
            if n == 3 and (a[0] != a[1] and a[0] != a[2] and a[1] != a[2]) and (cor11[0] != cor11[1]) and (
                    cor11[0] != cor11[2]) and (cor11[1] != cor11[2]) and (cor22[0] != cor22[1]) and (
                    cor22[0] != cor22[2]) and (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 (cor11[0] != cor11[1]) and (cor11[0] != cor11[2]) and (cor11[0] != cor11[3]) and (
                    cor11[1] != cor11[2]) and (cor11[1] != cor11[3]) and (cor11[2] != cor11[3]) and (
                    cor22[0] != cor22[1]) and (cor22[0] != cor22[2]) and (cor22[0] != cor22[3]) and (
                    cor22[1] != cor22[2]) and (cor22[1] != cor22[3]) and (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

LSM()函数实现最小二乘法估计，用于在特征匹配点基础上求解最小误差的仿射变换参数。

# LSM函数是为了使用最小二乘法来计算两组点之间的仿射变换参数。
import numpy as np
import math


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

该函数基于最大索引图，统计局部区域内方向响应直方图构建特征向量：

import numpy as np
import math
from einops import repeat

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.float64)
    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

为提升旋转不变性，针对输入图像进行角度迭代：

import sys
import os
from phasepack import phasecong
from FSC import fsc
from RIFT import rift
from RIFT_Descriptor_No_Rotation_Invariance import RIFT_descriptor_no_rotation_invariance
import numpy as np
import cv2

sys.path.append(os.path.join(os.path.dirname(__file__), "."))
sys.path.append(os.path.join(os.path.dirname(__file__), ".."))

step = 10
device = "cuda:0"


def rotation_optimization(img1, img2, best_mean_error, best_inliersIndex, best_match_point1, best_match_point2):
    best_rot_img = img2
    # 计算图像中心点，每次旋转angle度
    for angle in range(step, 360, step):
        print("=================================")
        print(f"旋转角度为{angle}")
        height, width = img1.shape[:2]
        center = (width // 2, height // 2)

        # 旋转源图像
        rotation_matrix = cv2.getRotationMatrix2D(center, angle, 1.0)
        rotated_img2 = cv2.warpAffine(img2, rotation_matrix, (width, height))

        # RIFT
        current_mean_error, inliersIndex, match_point1, match_point2 = rift(img1, rotated_img2)
        print(f"current_cost为{current_mean_error}")

        # 判断
        if current_mean_error < best_mean_error:
            best_mean_error = current_mean_error
            best_inliersIndex = inliersIndex
            best_match_point1 = match_point1
            best_match_point2 = match_point2
            best_rot_img = rotated_img2

    return best_inliersIndex,best_match_point1, best_match_point2, best_rot_img

三、总结

RIFT算法突破了传统图像匹配方法在多模态场景下的限制，提供了对辐射不一致性和局部结构信息提取的强鲁棒性。结合本文提供的Python实现代码，可以清晰了解RIFT从图像预处理到变换估计的全过程。

参考文献：J. Li, Q. Hu and M. Ai, “RIFT: Multi-Modal Image Matching Based on Radiation-Variation Insensitive Feature Transform,” in IEEE Transactions on Image Processing, vol. 29, pp. 3296-3310, 2020