Closest pair of a set of points in the plane
% closest(Points, Point1, Point2, DistSqu) is true if Point1 and Point2 are two
% points which are nearest together in the list of Points, and DistSqu is the
% square of the distance between the two points.
% e.g. closest([p(2,3),p(5,4),p(9,6),p(4,7),p(8,1),p(7,2)], P, Q, D)
% gives P=p(7,2), Q=p(8,1), D=2
closest(Points, Point1, Point2, DistSqu):-
Points=[_,_|_],
merge_sort(Points, 0, SortX), % Sorted according to the x-coordinate
length(SortX, Length),
closest_1(Length, SortX, [], t(Point1,Point2,DistSqu)).
% Note the similarity with merge_sort_1/5 below.
closest_1(1, [P|Rest], Rest, t(P,P,2147483647)):-!.
closest_1(2, [P,Q|Rest], Rest, t(P,Q,DistSqu)):-!,
distance_squared(P, Q, DistSqu).
closest_1(N, Left, Rest, Closest):-
N1 is N // 2,
N2 is N - N1,
closest_1(N1, Left, Right, ClosestL), % The closest pair in Left
closest_1(N2, Right, Rest, ClosestR), % The closest pair in Right
closer(ClosestL, ClosestR, ClosestLR), % The closest pair in Left or Right
Right=[p(Xmedian,_)|_], % The median of the N points
ClosestLR=t(_,_,DistSqu),
D=sqrt(DistSqu), % The distance between closest pair
Xmin is Xmedian - D,
left_strip(N1, Left, Xmin, StripL0), % Points > Median-D [and < Median]
Xmax is Xmedian + D,
right_strip(N2, Right, Xmax, StripR0), % Points [>= Median and] < Median+D
merge_sort(StripL0, 1, StripL), % Sorted according to the y-coordinate
merge_sort(StripR0, 1, StripR), % Sorted according to the y-coordinate
join(StripL, StripR, D, ClosestLR, Closest).
closer(ClosestA, t(_,_,E), ClosestA):-
ClosestA=t(_,_,D),
D < E, !.
closer(_, ClosestB, ClosestB).
% left_strip(N, Points, X, Left) is true if Left contains the elements amongst
% the first N elements of Points having an x-coordinate greater than X.
left_strip(0, _, _, []):-!.
left_strip(N, [Point|Points], X, [Point|Left]):-
Point=p(X1,_),
X1 > X, !,
N1 is N - 1,
left_strip(N1, Points, X, Left).
left_strip(N, [_|Points], X, Left):-
N1 is N - 1,
left_strip(N1, Points, X, Left).
% right_strip(N, Points, X, Right) is true if Right contains the elements
% amongst the first N elements of Points having an x-coordinate less than X.
right_strip(0, _, _, []):-!.
right_strip(N, [Point|Points], X, [Point|Right]):-
Point=p(X1,_),
X1 < X, !,
N1 is N - 1,
right_strip(N1, Points, X, Right).
right_strip(N, [_|Points], X, Right):-
N1 is N - 1,
right_strip(N1, Points, X, Right).
% Finds the closest pair, if any, with one point in Ls and the other in Rs.
join([], _, _, Soln, Soln):-!. % No more Ls; succeed
join(_, [], _, Soln, Soln):-!. % No more Rs; succeed
join([p(_,Yl)|Ls], Rs, D, Soln0, Soln):-
Rs=[p(_,Yr)|_],
Yl =< Yr - D, !, % L is below R-D
join(Ls, Rs, D, Soln0, Soln). % Try next L with same R
join(Ls, [p(_,Yr)|Rs], D, Soln0, Soln):-
Ls=[p(_,Yl)|_],
Yl >= Yr + D, !, % L is above R+D
join(Ls, Rs, D, Soln0, Soln). % Try same L with next R
join([L|Ls], Rs, D, Soln0, Soln):- % L is above R-D and below R+D
join_1(Rs, L, D, Soln0, Soln1), % Find improvements with the current L
join(Ls, Rs, D, Soln1, Soln). % Try next L with same R
% Finds the closest pair, if any, with one point being L and the other in Rs.
join_1([], _, _, Soln, Soln). % No more Rs; succeed
join_1([p(_,Yr)|_], p(_,Yl), D, Soln, Soln):-
Yl >= Yr + D, !. % L is above R+D; succeed
join_1([R|Rs], L, D, t(_,_,DistSqu0), Soln):- % L is below R+D
distance_squared(L, R, DistSqu1),
DistSqu1 < DistSqu0, !, % L and R form a closer pair
join_1(Rs, L, D, t(L,R,DistSqu1), Soln).
join_1([_|Rs], L, D, Soln0, Soln):- % L is below R+D
join_1(Rs, L, D, Soln0, Soln). % L and R are not a closer pair
% merge_sort(List, Axis, SortedList) is true if SortedList is the result of
% sorting the List of 2-dimensional points along the given Axis.
% e.g. merge_sort([p(3,1),p(4,1),p(5,9),p(2,6),p(5,3)], 0, X) gives
% X=[p(2,6),p(3,1),p(4,1),p(5,9),p(5,3)]
% merge_sort([p(3,1),p(4,1),p(5,9),p(2,6),p(5,3)], 1, X) gives
% X=[p(3,1),p(4,1),p(5,3),p(2,6),p(5,9)]
merge_sort(List, Axis, SortedList):-
length(List, Length),
merge_sort_1(Length, Axis, List, [], SortedList).
merge_sort_1(0, _, Rest, Rest, []):-!.
merge_sort_1(1, _, [A|Rest], Rest, [A]):-!.
merge_sort_1(N, Axis, List, Rest, Sorted):-
N1 is N // 2,
N2 is N - N1,
merge_sort_1(N1, Axis, List, TempList, SortedLeft),
merge_sort_1(N2, Axis, TempList, Rest, SortedRight),
ordered_merge(SortedLeft, SortedRight, Axis, Sorted).
ordered_merge([], Ys, _, Ys).
ordered_merge([X|Xs], [Y|Ys], Axis, [Y|Zs]):-
lt(Axis, Y, X), !,
ordered_merge([X|Xs], Ys, Axis, Zs).
ordered_merge([X|Xs], Ys, Axis, [X|Zs]):-
ordered_merge(Xs, Ys, Axis, Zs).
% Comparisons on a given Axis of two 2-dimensional points.
lt(0, p(P,_), p(R,_)):-P < R, !.
lt(1, p(_,Q), p(_,S)):-Q < S.
distance_squared(p(X1,Y1), p(X2,Y2), DistanceSquared):-
X1_X2 is X1 - X2,
Y1_Y2 is Y1 - Y2,
DistanceSquared is X1_X2*X1_X2 + Y1_Y2*Y1_Y2.
% naive_closest(Points, Point1, Point2, DistSqu) is true if Point1 and Point2
% are two points which are nearest together in the list of Points, and
% DistSqu is the square of the distance between the two points.
% e.g. naive_closest([p(2,3),p(5,4),p(9,6),p(4,7),p(8,1),p(7,2)], P, Q, D)
% gives P=p(8,1), Q=p(7,2), D=2
naive_closest(Points, Point1, Point2, DistSqu):-
naive_closest_1(Points, t(void,void,2147483647), t(Point1,Point2,DistSqu)).
naive_closest_1([], Closest, Closest).
naive_closest_1([Point|Points], Closest0, Closest):-
naive_closest_2(Points, Point, Closest0, Closest1),
naive_closest_1(Points, Closest1, Closest).
naive_closest_2([], _, Closest, Closest).
naive_closest_2([Point2|Points], Point1, t(_,_,DistSqu0), Closest):-
distance_squared(Point1, Point2, DistSqu1),
DistSqu1 =< DistSqu0, !,
naive_closest_2(Points, Point1, t(Point1,Point2,DistSqu1), Closest).
naive_closest_2([_|Points], Point1, Closest0, Closest):-
naive_closest_2(Points, Point1, Closest0, Closest).
% Generates N random points with X and Y values between Min and Max.
% e.g. rand_points(100, -999, 999, Ps), closest(Ps, P, Q, D)
% e.g. rand_points(100, -999, 999, Ps), naive_closest(Ps, P, Q, D)
rand_points(0, _, _, []):-!.
rand_points(N, Min, Max, [P|Ps]):-
N > 0,
rand_point(Min, Max, P),
N1 is N - 1,
rand_points(N1, Min, Max, Ps).
% Generates a random point with X and Y values between Min and Max.
rand_point(Min, Max, p(X,Y)):-
Min < Max,
rand_int(Min, Max, X),
rand_int(Min, Max, Y).
% rand_int(I, J, K) is true if K is a pseudo-random integer in the range I to J
% inclusive.
rand_int(I, J, K):-K is int(rand(J - I + 1)) + I.
% length(Xs, L) is true if L is the number of elements in the list Xs.
%length(Xs, L):-length_1(Xs, 0, L).
%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).
LPA Index
Home Page