Minimum Weight Triangulation (Dynamic Programming)

/* mwt(+Polygon, -Indices, -Coords, -Length) is true if Indices are the      */
/*   indices, and Coords are the coordinates of the endpoints of the         */
/*   diagonals such that the triangulation of the convex Polygon is optimal  */
/*   in the sense that the sum of the lengths of the diagonals, Length, is   */
/*   minimized. The vertices of the polygon must be given in sequence.       */
/* e.g. mwt([p(0,0),p(2,0),p(6,1),p(5,6),p(2,5),p(0,3)], Is, Cs, L).         */
/* gives Is=[e(2,6),e(2,5),e(3,5)]                                           */
/*       Cs=[d(p(2,0),p(0,3)), d(p(2,0),p(2,5)), d(p(6,1),p(2,5))]           */
/*   and L=14.262405525                                                      */
mwt(Polygon, Indices, Coords, Length):-
  length(Polygon, N),
  distances(Polygon),              % Initialises the database w(I,J,Wij)
  mwt_1(2, N, [], Length, Kss),    % Uses the database
  endpoints([[]|Kss], N, Indices), % Add row 1 to Kss
  diagonals(Indices, Polygon, Coords).
% for i=2 to n
% e.g. Called from mwt([p(0,0),p(2,0),p(6,1),p(5,6),p(2,5),p(0,3)], Is, Cs).
% gives L=14.262405525
%   and Kss=[[],[2],[3,3],[2,3,4],[2,5,5,5]]
mwt_1(I, N, Cols, Length, []):-
  I > N,
  N1 is N - 1,
  nth_member(Cols, N1, [T|_]),
  w(1,N,W), !,
  Length is T - W.
mwt_1(I, N, Cols, Length, [Ks|Kss]):-
  J is I - 2,
  mwt_2(Cols, [0], I, J, Cols1, [], Ks),
  I1 is I + 1,
  mwt_1(I1, N, [[0]|Cols1], Length, Kss).

% for j=i-1 downto 1
mwt_2([], _, _, 0, [], Kss, Kss).
mwt_2([ColJ|Cols0], RowI, I, J, [[Tij|ColJ]|Cols], Kss0, Kss):-
  K is I - 1,
  w(J,I,Wij), !,
  next_element(RowI, ColJ, K, Wij, 1.0E308, Tij, 0, Kij),
  % Tij is the length of the diagonals plus the length of the edge (i,j), of
  %   the minimum weight triangulation of the polygon with vertices j,j+1,...,i
  % Kij is the value of k such that the polygons j,j+1,...,k and k,k+1,...,i
  %   and the triangle j,k,i constitute the minimum weight triangulation of
  %   the polygon j,j+1,...,i
  J1 is J - 1,
  append(RowI, [Tij], RowI1),
  mwt_2(Cols0, RowI1, I, J1, Cols, [Kij|Kss0], Kss).

% for k=i-1 downto j+1
% Given elements T[I,J+1]...T[I,I-1], elements T[J+1,J]...T[I-1,J], 
%   K (=I-1), Wij and "infinity", calculates T[I,J] and K[I,J]
next_element([], [], _, _, Tij, Tij, Kij, Kij).
next_element([Tik|RowI], [Tkj|ColJ], K, Wij, Tij0, Tij, _, Kij):-
  Tij1 is Tik + Tkj + Wij,
  Tij1 < Tij0, !,
  K1 is K - 1,
  next_element(RowI, ColJ, K1, Wij, Tij1, Tij, K, Kij).
next_element([_|RowI], [_|ColJ], K, Wij, Tij0, Tij, Kij0, Kij):-
  K1 is K - 1,
  next_element(RowI, ColJ, K1, Wij, Tij0, Tij, Kij0, Kij).

% e.g. endpoints([[],[],[2],[3,3],[2,3,4],[2,5,5,5]], 6, Vss)
% gives Vss=[e(2,6),e(2,5),e(3,5)]
endpoints(Kss, N, Vss):-
  nth_member(Kss, N, [K|_]), !, % T[N,1] was derived from T[N,K] and T[K,1]
  endpoints_1([e(N,K),e(K,1)], Kss, Vss).

endpoints_1([], _, []).
endpoints_1([e(I,J)|IJs], Kss, Vss):-
  J is I - 1, !,                 % T[I,I-1] doesn't correspond to a diagonal
  endpoints_1(IJs, Kss, Vss).
endpoints_1([e(I,J)|IJs], Kss, [e(J,I)|Vss]):- % But T[I,J] does
  nth_member(Kss, I, Ks),
  nth_member(Ks, J, K), !,      % T[I,J] was derived from T[I,K] and T[K,J]
  endpoints_1([e(I,K),e(K,J)|IJs], Kss, Vss).

% Converts endpoints of the diagonals from indices to coordinates
% e.g. Vss=[e(2,6),e(2,5),e(3,5)],
%      Polygon=[p(0,0),p(2,0),p(6,1),p(5,6),p(2,5),p(0,3)],
%      diagonals(Vss, Polygon, Coords).
% gives Coords=[d(p(2,0),p(0,3)),d(p(2,0),p(2,5)),d(p(6,1),p(2,5))]
diagonals([], _, []).  
diagonals([e(I,J)|Vss], Polygon, [d(Pi,Pj)|Coords]):-
  nth_member(Polygon, I, Pi),
  nth_member(Polygon, J, Pj), !,
  diagonals(Vss, Polygon, Coords).

% Finds the distance D between each pair of points p(X1,Y1) and p(X2,Y2) in Ps
%   and asserts w(I,J,D), where I < J, and I and J are the indices of p(X1,Y1)
%   and p(X2,Y2) respectively in Ps
% e.g. distances([p(0,0),p(3,0),p(0,4)]).
% asserts w(1,2,3), w(1,3,4) and w(2,3,5)
  distances_1(Ps, 1).

distances_1([], _).
distances_1([P|Ps], I):-
  I1 is I + 1,
  distances_2(Ps, P, I, I1),
  distances_1(Ps, I1).

distances_2([], _, _, _).
distances_2([p(X1,Y1)|Ps], p(X2,Y2), I, J):-
  D is sqrt((X1-X2)*(X1-X2) + (Y1-Y2)*(Y1-Y2)),
  J1 is J + 1,
  distances_2(Ps, p(X2,Y2), I, J1).

/* append(Xs, Ys, Zs) is true if Zs is the result of appending the list Xs   */
/*   to the list Ys.                                                         */
%append([], Ys, Ys).
%append([X|Xs], Ys, [X|Zs]):-append(Xs, Ys, Zs).

/* 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(Xs, L0, L) is true if L is equal to L0 plus the number of        */
/*   elements in the list Xs.                                                */
%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

/* nth_member(+Xs, ?N, ?X) is true if X is the N-th (base 1) element of the  */
/*   list Xs. If N is uninstantiated, on backtracking all solutions will be  */
/*   found.                                                                  */
nth_member(Xs, N, X):-nth_member_1(Xs, X, 1, N).

nth_member_1([X|_], X, I, I).
nth_member_1([_|Xs], X, I0, I):-
  I1 is I0 + 1,
  nth_member_1(Xs, X, I1, I).

LPA Index     Home Page