cp_mismatchRemoval函数记录（使用修改后的RANSCA算法）

七仔的鸽子

已于 2023-04-16 18:20:01 修改

阅读量87

点赞数

分类专栏：图像匹配算法文章标签： matlab 算法

于 2023-04-13 17:27:37 首次发布

本文链接：https://blog.csdn.net/weixin_47000169/article/details/130135333

版权

图像匹配算法专栏收录该内容

10 篇文章 1 订阅

订阅专栏

输入：regis_points11,regis_points22,I1,I2,maxErr

输出：

function correctIndex = cp_mismatchRemoval(p1,p2,edge1,edge2,maxErr)
% mismatch removal algorithm
% Input: cp_mismatchRemoval(p1,p2,edge1,edge2,maxErr,iteration)
% p1: coarse match of image 1
% p2: coarse match of image 2
if nargin == 5
    iteration = 1;
    zoomflag = 1;
end
if length(p1) == 3
    correctIndex = [1 2 3];
    regis2sub = p2;
    return;
end
if length(p1) == 2
    correctIndex = [1 2];
    return;
end
correctIndex = [];
minArea1 = size(edge1,1) * size(edge1,2) / 80;
minArea2 = size(edge2,1) * size(edge2,2) / 80;

len = size(p1,1);
eor = [];
%% RANSAC algorithm
% if  ~zoomflag
%     for i = 1:len-2
%         for j = i+1:len-1
%             for k = j+1:len
%                 pp1 = [p1(i,1:2);p1(j,1:2);p1(k,1:2)];
%                 pp2 = [p2(i,1:2);p2(j,1:2);p2(k,1:2)];
%                 A1 = p1(i,:)-p1(j,:);
%                 B1 = p1(i,:)-p1(k,:);
%                 A2 = p2(i,:)-p2(j,:);
%                 B2 = p2(i,:)-p2(k,:);
%                 if 4*abs( A1(1)*B1(2)-A1(2)*B1(1) ) < minArea1 || 4*abs( A2(1)*B2(2)-A2(2)*B2(1) ) < minArea2
%                     continue;
%                 end
%                 side1A = norm(A1); side1B = norm(B1);
%                 side2A = norm(A2); side2B = norm(B2);
%                 A2A = side1A / side1B;
%                 B2B = side2A / side2B;
%                 if  A2A/B2B < 0.92 || A2A/B2B > 1.08
%                     continue;
%                 end
%                 ERR = getEdgeOverlappingRatio(pp1,pp2,p1(:,1:2),p2(:,1:2),i,j,k,maxErr);
%                 eor = [eor;i j k ERR]; 
%             end
%         end
%     end
% else
    for n = 1:min(len*(len-1)*(len-2)/6,500) 
        ijk = randperm(len,3);
        ir = ijk(1);
        jr = ijk(2);
        kr = ijk(3);
        pp1 = [p1(ir,1:2); p1(jr,1:2); p1(kr,1:2)];
        pp2 = [p2(ir,1:2); p2(jr,1:2); p2(kr,1:2)];
        A1 = p1(ir,:)-p1(jr,:);
        B1 = p1(ir,:)-p1(kr,:);
        A2 = p2(ir,:)-p2(jr,:);
        B2 = p2(ir,:)-p2(kr,:);
        if abs( A1(1)*B1(2)-A1(2)*B1(1) ) < minArea1 || abs( A2(1)*B2(2)-A2(2)*B2(1) ) <minArea2
            continue;
        end
        ransacERR = getEdgeOverlappingRatio(pp1,pp2,p1(:,1:2),p2(:,1:2),ir,jr,kr,maxErr);
        eor = [eor; ir jr kr ransacERR];
    end
% end
%% 
if size(eor,1) < 3
    error('No sufficent matches! ( less than 3 matches obtained )');
end
[~,ind] = max(eor(:,5));
base1 = p1(eor(ind,1:3),1:2);
base2 = p2(eor(ind,1:3),1:2);
affmat0 = fitgeotrans(base1,base2,'affine');
correctIndex = eor(ind,1:3)';
for i = 1:len
    if i == eor(ind,1) | i == eor(ind,2) | i == eor(ind,3)
        continue;
    end
    pp1_aff = [p1(i,1:2) 1] * affmat0.T;
    pp2_to_pp1aff = abs(p2(i,1:2) - pp1_aff(1:2));
    if (pp2_to_pp1aff(1) < maxErr / 1.5 & pp2_to_pp1aff(2) < 1.5 * maxErr) | ...
           (pp2_to_pp1aff(1) < 1.5 * maxErr & pp2_to_pp1aff(2) < maxErr / 1.5)
        correctIndex = [correctIndex; i];
    end
end
end
%% RANSAC algorithm
function RMSE = getEdgeOverlappingRatio(pp1, pp2, p1, p2, a, b, c, maxErr)
affmat = fitgeotrans(pp1,pp2,'affine');
p1_aff = [p1 ones(length(p1),1)] * affmat.T;

p2_to_p1aff = (p2 - p1_aff(:,1:2)).^2; % Distance between origin points and affined points
RMSE(1,1) = sqrt(sum(p2_to_p1aff(:)) / length(p2));

pixelDistance = p2_to_p1aff(:,1) + p2_to_p1aff(:,2);
pixelDistance([a b c]) = 2*maxErr^2; % ignore origin three points
[u,~] = find(pixelDistance <  maxErr^2); % mount of internal points of fitting result
[minerr,ind0] = min(pixelDistance );
RMSE(1,2:4) = [length(u) minerr ind0];
end

第一部分：初始化变量

除以80的目的是为了计算图像边缘的最小区域值。具体而言，该代码首先获取了第一幅图像和第二幅图像边缘矩阵的尺寸（即行数和列数），然后将两个尺寸相乘并除以80，计算出每个图像边缘所占用的最小区域。

这一操作的目的是为了避免将太小的边缘区域误判为特征点。由于较小的边缘区域可能只代表图像中某些微小的纹理或噪声，而非真正的特征，因此需要设置一个最小区域阈值来排除这些无意义的边缘区域。

具体的数值80是通过经验得出的，不同的实验数据和实际应用场景下可能会有不同的取值，需要根据具体情况进行调整。

correctIndex = [];
minArea1 = size(edge1,1) * size(edge1,2) / 80;
minArea2 = size(edge2,1) * size(edge2,2) / 80;

len = size(p1,1);
eor = [];

第二部分：

每一幅图片选出随机的3个点，计算其由对角线上的点所构成的矩形的面积，并将它们与最小面积进行对比，筛选出大于最小面积的点。

将符合要求的点带入到getEdgeOverlappingRatio函数进行计算。

for n = 1:min(len*(len-1)*(len-2)/6,500) 
    ijk = randperm(len,3);
    ir = ijk(1);
    jr = ijk(2);
    kr = ijk(3);
    pp1 = [p1(ir,1:2); p1(jr,1:2); p1(kr,1:2)];
    pp2 = [p2(ir,1:2); p2(jr,1:2); p2(kr,1:2)];
    A1 = p1(ir,:)-p1(jr,:);
    B1 = p1(ir,:)-p1(kr,:);
    A2 = p2(ir,:)-p2(jr,:);
    B2 = p2(ir,:)-p2(kr,:);
    if abs( A1(1)*B1(2)-A1(2)*B1(1) ) < minArea1 || abs( A2(1)*B2(2)-A2(2)*B2(1) ) <minArea2
        continue;
    end
    ransacERR = getEdgeOverlappingRatio(pp1,pp2,p1(:,1:2),p2(:,1:2),ir,jr,kr,maxErr);
    eor = [eor; ir jr kr ransacERR];
end

下面就是getEdgeOverlappingRatio函数。

%% RANSAC algorithm
function RMSE = getEdgeOverlappingRatio(pp1, pp2, p1, p2, a, b, c, maxErr)
affmat = fitgeotrans(pp1,pp2,'affine');
p1_aff = [p1 ones(length(p1),1)] * affmat.T;

p2_to_p1aff = (p2 - p1_aff(:,1:2)).^2; % Distance between origin points and affined points
RMSE(1,1) = sqrt(sum(p2_to_p1aff(:)) / length(p2));

pixelDistance = p2_to_p1aff(:,1) + p2_to_p1aff(:,2);
pixelDistance([a b c]) = 2*maxErr^2; % ignore origin three points
[u,~] = find(pixelDistance <  maxErr^2); % mount of internal points of fitting result
[minerr,ind0] = min(pixelDistance );
RMSE(1,2:4) = [length(u) minerr ind0];
end

其中：

tform = fitgeotrans(movingPoints,fixedPoints,transformationType)

movingPoints — 图像上想要移动的点的坐标，至少是两个double型2维点.
fixedPoints — 目标点，和上面同等规模
transformationType — 变换类型，包括如下几种：

transformationType	Description
‘Affine’	仿射变换
‘NonreflectiveSimilarity’	非反射相似变换
‘Projective’	投影变换
‘Similarity’	相似变换（即仿射变换中去除错切变换）

简单来说就是将pp1经过仿射变换移动到pp2处，返回变换矩阵。

在本例中，使用3个点进行仿射变换，得到变换矩阵之后，给总点矩阵乘以变换矩阵得到将p1全部的点变换到p2处对应的点。计算与原始p2处的点的欧式距离，保存到返回变量RMSE中。

由于 RANSAC 算法对一些不良的匹配点数据敏感，因此该函数将pixelDistance 中前三个特征点的距离值直接抛弃，避免它们对后续结果造成干扰。

筛选出距离<maxeErr^2的距离个数，最小的距离与最小距离的索引存入RMSE 。

第三部分：

if size(eor,1) < 3
    error('No sufficent matches! ( less than 3 matches obtained )');
end
[~,ind] = max(eor(:,5));
base1 = p1(eor(ind,1:3),1:2);
base2 = p2(eor(ind,1:3),1:2);
affmat0 = fitgeotrans(base1,base2,'affine');
correctIndex = eor(ind,1:3)';
for i = 1:len
    if i == eor(ind,1) | i == eor(ind,2) | i == eor(ind,3)
        continue;
    end
    pp1_aff = [p1(i,1:2) 1] * affmat0.T;
    pp2_to_pp1aff = abs(p2(i,1:2) - pp1_aff(1:2));
    if (pp2_to_pp1aff(1) < maxErr / 1.5 & pp2_to_pp1aff(2) < 1.5 * maxErr) | ...
           (pp2_to_pp1aff(1) < 1.5 * maxErr & pp2_to_pp1aff(2) < maxErr / 1.5)
        correctIndex = [correctIndex; i];
    end
end

选出上一步中处理得到的RMSE中距离<maxErr^2最多的值，按这个值选出eor中这三个随机数，再通过随机数选出p1与p2中的x，y值赋予base1，base2。

将base1通过仿射变换映射到base2处，得变换矩阵（此时的变换矩阵是相对于其他最精确的），之后再计算曼哈顿距离，筛选出符合范围的点的索引。