RIFT Multi-Modal Image Matching Based on Radiation-Variation Insensitive Feature Transform

RIFT:radiation-variation insensitive feature transform

introduction

3个贡献：
1.使用时相一致性而不是像素强度做特征检测，更鲁棒；只检测角点和边缘点
2.提出 maximum index map来做特征描述
3.支持旋转不变性

要实现通用的鲁棒图像匹配，需要解决三个难题:(1)几何变换不变性，辐射不变性（2）外点的去除（3）非刚性的配准问题
RIFT主要用来解决辐射不变性，尤其是NRD(nonlinear radiation distortiion).
目前，如果不知道图像的几何信息，没有任何一个方法可以同时用于optical-optical,sar-optial,day-night等多模态的配准。

在PC map的角点和边缘点检测特征，这些点的鲁棒性较好
使用log-Gabor构建maximum index map.
构建multiple MIM来解决旋转不变性。
HOPC的不足：（1）需要提前知道图像的地理信息，利用地理信息做粗匹配。但地理信息有可能本身差异很大（2）HOPC依赖模板匹配来寻找对应点，对于旋转和尺度变换不适用（3）HOPC使用Harris寻找关键点，关键点可能找的不好
RIFT不依赖地理信息，具有旋转不变性

method

在两张图上检测关键点；关键点的特征提取；通过特征匹配来匹配关键点；外点过滤；homograpy计算
所以，关键点的鲁棒性很重要（即同一个关键点可以在不同的图像中都被检测到）。

1.角点和边缘点检测：
输入一对图像（光学图和雷达图）（a），直接使用FAST算法检测两张图像的角点和边缘点(d)，效果不好（两张图的角点并不对应，雷达图中很多角点和边缘都没有检测到）。因此，我们利用PC来检测角点。其过程为：
首先，利用小波变换得到光学图和雷达图的PC图。小波变换可以获得 freq-time图，PC图即计算各个波间相位的一致性。

然后，根据PC图计算moment,得到minimum moment和maximum moment.。minimum moment大的地方对应角点（minimum moment经过极大值检测和非极大值抑制后得到的点，记作角点，图b,e）。然后用FAST算法对maximum moment做边缘检测，得到边缘点（图ef）.

对应的代码为：

[des_m1,des_m2] = RIFT_no_rotation_invariance(im1,im2,4,6,96);#在两张图像中都就检测关键点，并且将这些关键点的特征提取出来[m1,~,~,~,~,eo1,~] = phasecong3(im1,s,o,3,'mult',1.6,'sigmaOnf',0.75,'g', 3, 'k',1);#m1是maximum moment[m2,~,~,~,~,eo2,~] = phasecong3(im2,s,o,3,'mult',1.6,'sigmaOnf',0.75,'g', 3, 'k',1);#m2是另一张图的maximum momenta=max(m1(:)); b=min(m1(:)); m1=(m1-b)/(a-b);#根据最大最小值做归一化a=max(m2(:)); b=min(m2(:)); m2=(m2-b)/(a-b);% FAST detector on the maximum moment maps to extract edge feature points. 使用FAST算法在m1,m2图中检测边缘点，找出5000个点m1_points = detectFASTFeatures(m1,'MinContrast',0.05);m2_points = detectFASTFeatures(m2,'MinContrast',0.05);m1_points=m1_points.selectStrongest(5000);   %number of keypoints can be set by usersm2_points=m2_points.selectStrongest(5000);

2.特征提取
先构造MIM图，然后在MIM图中，按照SIFT的特征提取算法来提取特征。
先当与先将两幅图转为同一种图，然后在该图中提取特征做配准。
(1) MIM的构造：

代码：

RIFT_no_rotation_invariance(im1,im2,4,6,96);
RIFT_descriptor_no_rotation_invariance(im1, m1_points.Location,eo1, patch_size, s,o）
其中s=4,o=6，eo1是一个4x6的cell,每个元素是HxW(与图像尺寸一致)的数组，表示该scale,该orientation下对应的小波变换与原图做卷积的结果。

关键：
CS(:,:,j)=CS(:,:,j)+abs(eo{i,j});#将各个scale的变换结果的幅度相加（一张图的一个尺度一个角度的小波变换结果是复数平面，abs取复数的幅度）

[~, MIM] = max(CS,[],3); % MIM maximum index map 在6个角度中，查看哪个角度的幅度最大；MIM为幅度最大的角度的索引

function des = RIFT_descriptor_no_rotation_invariance(im, kps,eo, patch_size, s,o)KPS=kps'; %keypoints
[yim,xim,~] = size(im);CS = zeros(yim, xim, o); %convolution sequence
#将各个scale的变换结果图相加，其中scale对应频率。
for j=1:ofor i=1:sCS(:,:,j)=CS(:,:,j)+abs(eo{i,j});end
end
#CS的shape为[H,W,6],6表示6个方向，max(CS,[],3)表示返回CS中第3个维度的最大值。所以MIM的shape为H，W
[~, MIM] = max(CS,[],3); % MIM maximum index map

(2)采用类似于SIFT的方法，使用MIM图构造特征
SIFT构造特征的方法：先把16x16的区域分为16个4x4的小区域，每个小区域中汇总16个像素的梯度方向（每个像素的梯度强度根据其方向分配到8个方向中）作为该区域的特征向量。

RIFT的patchsize是96x96，分为了6x6个小区域（每个小区域的尺寸为16x16），然后在每个小区域内统计值的直方图（每个像素的的值按照其方向分配到6个方向中，这里没有计算梯度的步骤，因为此时的值本身有梯度和方向的信息），每个小区域的特征向量长度为6.所以，一个patch的特征向量长度为6x6x6.
代码为：


des = zeros(36*o, size(KPS,2)); %descriptor (size: 6×6×o) 每个patch的特征向量的长度为36*6
kps_to_ignore = zeros(1,size(KPS,2)); #size(KPS,2)表示KPS的个数，即关键点的个数for k = 1: size(KPS,2)x = round(KPS(1, k));y = round(KPS(2, k));#x1,y1,x2,y2表示以该特征点为中心的patch的左上，右下角坐标x1 = max(1,x-floor(patch_size/2));y1 = max(1,y-floor(patch_size/2));x2 = min(x+floor(patch_size/2),size(im,2));y2 = min(y+floor(patch_size/2),size(im,1));    patch = MIM(y1:y2,x1:x2); %local MIM patch for feature description#将patch分为6*6的区域，每个区域单独提取特征ns=6;RIFT_des = zeros(ns,ns,o);  %descriptor vector% histogram vectorsfor j = 1:nsfor i = 1:nsclip = patch(round((j-1)*ys/ns+1):round(j*ys/ns),round((i-1)*xs/ns+1):round(i*xs/ns));RIFT_des(j,i,:) = permute(hist(clip(:), 1:o), [1 3 2]); #hist是对每个小区域提取特征，hist是统计值的分布，permute是改变array的shapeendendRIFT_des=RIFT_des(:);#特征做归一化if norm(RIFT_des) ~= 0RIFT_des = RIFT_des /norm(RIFT_des);enddes(:,k)=RIFT_des;
end
des = struct('kps', KPS(:,kps_to_ignore ==0)', 'des', des(:,kps_to_ignore==0)');

(2) 特征点匹配，与sift的方法一致（特征距离最小的认为匹配）

#matlab的matchFeatures中的默认匹配方法是Exhaustive,会一一计算输入特征组1和输入特征组2的两两指尖的特征距离，默认特征间的距离计算方式为SSD。
[indexPairs,matchmetric] = matchFeatures(des_m1.des,des_m2.des,'MaxRatio',1,'MatchThreshold', 100);
matchedPoints1 = des_m1.kps(indexPairs(:, 1), :);
matchedPoints2 = des_m2.kps(indexPairs(:, 2), :);
[matchedPoints2,IA]=unique(matchedPoints2,'rows'); #unique类似python的set，去掉重复的。rows表示每一行作为一个元素；返回结果为去重复的matchedPoints2和对应的索引
matchedPoints1=matchedPoints1(IA,:);

(3) 外点去除
先计算出Homography,然后根据该hoography判断已经配准的点哪些是inliners,从而在显示配准结果时将inliears和outliers分开显示

H=FCS(matchedPoints1,matchedPoints2,'affine',2);
Y_=H*[matchedPoints1';ones(1,size(matchedPoints1,1))];
Y_(1,:)=Y_(1,:)./Y_(3,:);
Y_(2,:)=Y_(2,:)./Y_(3,:);
E=sqrt(sum((Y_(1:2,:)-matchedPoints2').^2));
inliersIndex=E<3;
cleanedPoints1 = matchedPoints1(inliersIndex, :);
cleanedPoints2 = matchedPoints2(inliersIndex, :);

经过特征点匹配后得到的配准关系（cor1,cor2）有很多的outliears，不能直接计算homography,需要先判定哪些是outliers，哪些是inliears。这里采用的方法是：
（1）每次选择3个点，计算affine homography的6个参数
(2)根据计算得到的affine homograpy,判断内点的个数。
（3）不断迭代过程1和过程2，内点个数最多时的homography就是最优解
注意：这里的outliears和inliears判断时，文章中说采用NBSC方法，代码中说采用FSC方法；但实际的代码却与NBSC,FSC方法都不同。代码的思想与FSC,NBSC相似，都是遍历所有可能的affine model,内点最多的affine model为最优解。

function [solution,rmse,cor1_new,cor2_new]=FSC(cor1,cor2,change_form,error_t)
[M,N]=size(cor1); %M=1389,N=2,其中N=2表示坐标n=3; #表示计算affine homograpy只需要3对配准点iterations=10000; #表示最多迭代计算10000次
most_consensus_number=0; #初始化most_consensus_number为0，most_consensus_number表示inlieras个数的最大值
cor1_new=zeros(M,N); %（1389，2）
cor2_new=zeros(M,N);%（1389，2）for i=1:1:iterationsa=floor(1+(M-1)*rand(1,n)); #rand(1,n)随机得到0-1之间的数字n次；所以floor(1+(M-1)*rand(1,n))随机获得n次（1，M-1）之间的数#随机挑选了n个匹配点（n=3）cor11=cor1(a,1:2); cor22=cor2(a,1:2);#根据挑选的点计算affine的参数           [parameters,~]=LSM(cor11,cor22,change_form);solution=[parameters(1),parameters(2),parameters(5);parameters(3),parameters(4),parameters(6);parameters(7),parameters(8),1];#根据affine model来计算当前配准点的误差diff_match2_xymatch1_xy=cor1(:,1:2)';match1_xy=[match1_xy;ones(1,M)];t_match1_xy=solution*match1_xy;match2_xy=cor2(:,1:2)';match2_xy=[match2_xy;ones(1,M)];diff_match2_xy=t_match1_xy-match2_xy;diff_match2_xy=sqrt(sum(diff_match2_xy.^2));#误差小于阈值的记inliears,inliears的个数记作consensus_numindex_in=find(diff_match2_xy<error_t);consensus_num=size(index_in,2);#私用内点个数最多的情况，来计算affine homographyif(consensus_num>most_consensus_number)most_consensus_number=consensus_num;cor1_new=cor1(index_in,:);cor2_new=cor2(index_in,:);end
end
[parameters,rmse]=LSM(cor1_new(:,1:2),cor2_new(:,1:2),change_form);
solution=[parameters(1),parameters(2),parameters(5);parameters(3),parameters(4),parameters(6);parameters(7),parameters(8),1];
end

注：
NBSC算法过程（对应文章ROBUST FEATURE MATCHING FOR GEOSPATIAL IMAGES VIA AN AFFINE-INVARIANT COORDINATE）：
NBCS也是通过循环来遍历一个个可能的homography,然后inliears个数最多的homograpy对应最终的解。

其中NBCs表示在affine变换时，4个点构成的4个三角形，其面积比保持不变。如果不满足这个性质，则说明这4点的对应关系有误。

FSC算法过程：

result

关键点检测
检测到的关键点的可重复性：即同名点在两幅图中都作为关键点被检测到的能力。RIFT方法的效果最好