Matlab + Ransac 求两个对应点集的变换矩阵

Ransac的算法流程如下：

% ---- run RANSAC  ----
% define RANSAC parameters
threshDistBetPoints = 0.1;
numOfIteration = 50;
numOfRandomPoints = 3;
inlierRatio = 0.6;

% output parameters of RANSAC
best_transform = [];
best_num_inliers = 0;

% Plot the data points
for i=1:numOfIteration
%    1) select 3 random points
        idx = randperm(size(bestMatches,2),numOfRandomPoints); 
        sample = bestMatches(:,idx); 

%    2) compute affine transform T'

%    3) take all 500 matching points
        % map all Q(right image) to P'(left image) (transformation)
        % Compute distances between P' and P point (Left Image) 
        % compute no. of inliers
%    4) store the currently best affine transform
        % output of first iteration - T' and no. of inliers     
        if inlierNum >= best_num_inliers 
             best_num_inliers =  inlierNum
        if end
for end
% End of algorithm best transform = T' with largest numbers of inliers
% stitch function - already implemented

注意，Ransac输入的是两个已经匹配好的，对应起来的点，Ransac是在一堆大多数匹配好的点对寻找变换矩阵

下面是matlab代码：

% 输入长度相同的两个二维列向量矩阵
% 输出拟合的point1到point2的单应矩阵HH
% dd_num迭代次数
% theah是阈值，小于阈值认为是一个点
% Res_co是point1经过矩阵HH映射到point2的点
function [HH, Res_co] = findhomobyransac(point1,point2,dd_num,theah)
Res_co = [];
pot1_num = length(point1);
pot2_num = length(point2);

num = 0;
i=1;
best_num = 0;
if pot1_num == pot2_num && pot2_num > 3
    while i < dd_num
        temp = randperm(pot1_num,4);
        point11 = [point1(temp(1),:); point1(temp(2),:);point1(temp(3),:);point1(temp(4),:)];
        point22 = [point2(temp(1),:); point2(temp(2),:);point2(temp(3),:);point2(temp(4),:)];

        H = findhomeM(point11, point22);
        
        %如果求不出相应的矩阵，则直接结束此次循环
        if H(3,3) == 0
            continue
        end
        
        %返回映射过去的点的坐标
        Res_co = [];
        pre_num = 0;
        for j = 1: length(point1)
            % 求映射过去的点的坐标
            Coor = [(H(1,1)*point1(j,1) + H(1,2)*point1(j,2) + H(1,3))/(H(3,1)*point1(j,1) + H(3,2)*point1(j,2) + H(3,3))
                (H(2,1)*point1(j,1) + H(2,2)*point1(j,2) + H(2,3))/(H(3,1)*point1(j,1) + H(3,2)*point1(j,2) + H(3,3))
                1];
            Res_co = [Res_co; Coor'];
            % 遍历所有点，求在阈值内满足要求的点
            for jj = 1:length(point2)
            	d = sqrt((Coor(1)-point2(jj,1))^2 + (Coor(2)-point2(jj,2))^2);
                if d < theah
                    pre_num = pre_num + 1;
                    break
                end
            end
        end
        
        % 如果这次迭代结果比之前的都好，就把这次的结果赋值存下来
        if pre_num > best_num
            best_num = [pre_num,i];
            num = 0;
            HH = H;
        else
            num = num + 1;
        end
        
        % 迭代20次都有没有优化，就结束
        if num > 20
            break
        end
        i = i + 1;
%         best_num
    end
else
    HH = [0 0 0;0 0 0;0 0 0];
end

end

% 求解point到point1的变换矩阵，使用四个对应的点计算
function H = findhomeM(point, point1)

x1 = point(1,1);x2 = point(2,1);x3 = point(3,1);x4 = point(4,1);
y1 = point(1,2);y2 = point(2,2);y3 = point(3,2);y4 = point(4,2);

a1 = point1(1,1);a2 = point1(2,1);a3 = point1(3,1);a4 = point1(4,1);
b1 = point1(1,2);b2 = point1(2,2);b3 = point1(3,2);b4 = point1(4,2);

A = [x1, y1, 1, 0, 0, 0, -a1*x1, -a1*y1;
    0, 0, 0, x1, y1, 1, -b1*x1, -b1*y1;
    x2, y2, 1, 0, 0, 0, -a2*x2, -a2*y2;
    0, 0, 0, x2, y2, 1, -b2*x2, -b2*y2;
    x3, y3, 1, 0, 0, 0, -a3*x3, -a3*y3;
    0, 0, 0, x3, y3, 1, -b3*x3, -b3*y3;
    x4, y4, 1, 0, 0, 0, -a4*x4, -a4*y4;
    0, 0, 0, x4, y4, 1, -b4*x4, -b4*y4  
    ];
b = [a1; b1; a2; b2; a3; b3; a4; b4];
if rank(A) < 8
    H = [0 0 0;0 0 0;0 0 0];
else
    h = A \ b ;
    H = [h(1),h(2),h(3);h(4),h(5),h(6);h(7),h(8),1];
end
end

一般情况下总能找到正确的变换矩阵，但是本人觉得效果还是不如意，希望大家能够指出存在的问题并有所改进，有更好的方法记得分享给我372774377@qq.com