分治法和蛮力法MATLAB求最近点对

主程序 main.m

%%
clear;clc
n = 20;
% 随机生成20个点
A=rand(n,2)*10; 
% 将20个点按横坐标升序排列
A = sortrows(A,1);
% 蛮力法求随机点的最近点对
[mindist,x1,x2] = Bcloest(A,1,n);
[mindist1,y1,y2] = cloest(A,1,n);
% 使用红色的点标记随机点
x=A(:,1); 
y=A(:,2);
for i = 1:n
    plot(x,y,'r.'); hold on
    text(A(i,1),A(i,2),num2str(i));hold on
end
%使用绿色十字标记分治法的最近点对
plot(y1(1),y1(2),'g+');hold on
plot(y2(1),y2(2),'g+');hold on
% 使用蓝色的十字标记最近点对
plot(x1(1),x1(2),'b*');hold on
plot(x2(1),x2(2),'b*');hold on

分治求最近点对 cloest.m
function [d,x1,x2] = cloest(S,low,high)
P=[];
if(high - low == 1)
    [d,x1,x2] = deal(pdist2(S(low,:),S(high,:)),S(low,:),S(high,:));
elseif(high - low == 2)
    d1 = pdist2(S(low,:),S(low + 1,:));
    d2 = pdist2(S(low + 1,:),S(high,:));
    d3 = pdist2(S(low,:),S(high,:));
    if((d1 < d2) && (d1 < d3))
        [d,x1,x2] = deal(d1,S(low,:),S(low + 1,:));
    elseif(d2 < d3)
        [d,x1,x2] = deal(d2,S(low+1,:),S(high,:));
    else
        [d,x1,x2] = deal(d3,S(low,:),S(high,:));
    end
else
    mid = floor((low+high)/2);
    [d1,x11,x12] = cloest(S,low,mid);
    [d2,x21,x22] = cloest(S,mid+1,high);
    if(d1 <= d2)
        [d,x1,x2] = deal(d1,x11,x12);
    else
        [d,x1,x2] = deal(d2,x21,x22);
    end
    index = 0;
    for i = mid:-1:low
        if(S(mid,1) - S(i,1) >= d)
            break;
        end
        index = index + 1;
        P(index,:) = S(i,:);
     
    end
    for i = mid+1:high
        if(S(i,1) - S(mid,1) >= d)
            break;
        end
        index = index + 1;
        P(index,:) = S(i,:);
        
    end
     sortrows(P,2);
    for i = 1:index
        for j = i+1:index
            if(P(j,2) - P(i,2) >= d)
                break;
            else
                d3 = pdist2(P(i,:),P(j,:));
                if(d3 < d)
                    [d,x1,x2] = deal(d3,P(i,:),P(j,:));
                end
            end
        end
    end
    end

蛮力法求最近点对 Bcloest.m

function [d,x1,x2] = Bcloest(A,low,high)
min = 100;
for i = low:high
    for j = i + 1:high
        d1 = pdist2(A(i,:),A(j,:));
        if( d1 < min)
            [d,x1,x2] = deal(d1,A(i,:),A(j,:));
            min = d1;
        end
    end
end
end