my_ransac.m

This is a computer vision project implemented in matlab to remove radial distortion from the image

function [F,inlrs] = my_ransac(x,xp,dThresh)
if nargin<3
    dThresh = 1.5;
end
e = .60;                % outlier probability              
p = .99;
s = 8;
N = log(1 - p)/log(1 - (1 - e)^s);
samplesTaken = 0;
bestResErr = Inf;
maxInliers = 0;
while samplesTaken < N
    %select samples
    samples = randperm(size(x,1));
    samples = samples(1:15);
    samplesTaken = samplesTaken + 1;
   
    %compute model    
    xs = x(samples,1:2); xsp = xp(samples,1:2);
    %F = getF7pt( x(samples,1:2), xp(samples,1:2) );
    F = getF15pt2D( x(samples,1:2), xp(samples,1:2) );
    Fp = F(2:end, 2:end);
    
    six = F(2,1)/F(2,4);
    siy = F(1,2)/F(4,2);
    
    %determine inliers
    inliers = [];
    resErr = 0;
    for i=1: size(x)
        L1 = Fp*x(i,:)';
        %dist1 = abs(dot( xp', L1 ));
        Cx = [-L1(1,1)/(2*six*L1(3,1)), -L1(2,1)/(2*six*L1(3,1))];
        temp = L1(2,1)^2+L1(1,1)^2 - 4*six*L1(3,1);
        temp = temp/((six^2)*L1(3,1)^2);
        rx = 0.5*sqrt(temp);

        L2 = Fp' * xp(i,:)';
        Cy = [-L2(1,1)/(2*siy*L2(3,1)), -L2(2,1)/(2*siy*L2(3,1))];
        %dist2 = abs(dot( L2, x' ));
        temp = L2(2,1)^2+L2(1,1)^2 - 4*siy*L2(3,1);
        temp = temp/((siy^2)*L2(3,1)^2);
        ry = 0.5*sqrt(temp);

        % distance applies here
        err = (distan(x(i,1:2), Cx)-rx)^2+(distan(xp(i,1:2), Cy)-ry)^2;

        if err < dThresh
            inliers = [inliers,i];
        end
        resErr = resErr+err;
    end



    %inliers = find( dist1 < dThresh & dist2 < dThresh );
    inlierCount = size(inliers,2);            
    if inlierCount > 0
        resErr = resErr / inlierCount;
        %resErr = sampsonErrf( F, x(inliers,1:2), xp(inliers,1:2) );
    else
        resErr = Inf;
    end
    if inlierCount > maxInliers | ...
       (inlierCount == maxInliers & resErr < bestResErr)
        % keep best found so far
        maxInliers = inlierCount;
        bestResErr = resErr;
        bestF = F;
        bestInliers = inliers;

        % adaptively update N
        e = 1 - inlierCount / size( pmatches,1 );
        if e > 0
            N = log(1 - p)/log(1 - (1 - e)^s);
        else
            N = 1;
        end              
    end
end 

inlrs = bestInliers;
F = bestF;