输入: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处,得变换矩阵(此时的变换矩阵是相对于其他最精确的),之后再计算曼哈顿距离,筛选出符合范围的点的索引。