Zachary Danziger

Arthur,
If you are just looking for the standard Euclidean distance you can leave off the dfcn argument and it will be used by default.

Otherwise it expects a function handle. The 'ordinary' distance function is passed as:
dfcn = @(p,q) sqrt(sum( (p-q).^2 ));

But any function will work that operates on points p and q that respects their dimensionality. A Chebyshev-like measure could be used, for example:
dcheb = @(p,q) max(abs(p-q))

Lingji,

I cannot see a way to calculate the coupling distance during the recursive calls to c(i,j) without additional overhead. I feel that adding computational cost is not justified because, in general, the coupling sequence is not unique and therefore is not especially informative. The coupling sequence is just any allowable (i.e., follows the forward movement rules) sequence between points on P and Q that never has distance greater than cm. In the typical case there are a great many allowable sequences because interchanging many points with small distances does not affect cm if there is a much larger distance later in the sequence.

Here is my compromise: if the user requests the coupling sequence then I calculate the cm as usual, and at the end of the code I loop through the CA variable choosing one workable coupling sequence. I hope this pleases the fans.

Lingji,

I misunderstood your comment. Yes, the extra inputs to the c function are not needed. I have resubmitted the file with these removed for greater speed. Thank you for your suggestion.

Zach and Lingi; I would like to compare two drifter tracks (lat, long pairs). I tried using 'distance' from the Mapping Toolbox, but this dosn't work. Can you provide an example of dfcn function that does work? Thanks,

Zachary,

I hope you don't mind my posting this modified version; if you like it, feel free to use it to update your package or do whatever you see fit.

function [dFrechet, coupling] = fcn_discreteFrechetDistance(P, Q, dist_fcn)

% P and Q: two curves represented by matrices of size dim by number
% dist_fcn: optional distance function handle; default to Eucleandian

[n1, p] = size(P);
[n2, q] = size(Q);

if (n1 ~= n2)
error('P and Q mismatch');
end

switch nargin
case 2
dist_fcn = @(x, y) sqrt(sum((x - y).^2));
case 3
otherwise
error('wrong number of arguments');
end

% This is worked on by the inner function
CA = -ones(p, q);

% This is also worked on by the inner functions.
% Each cell stores four numbers [[u_last; v_last], [u_new; v_new]].
% The first column is an index into the cell array, and the second column
% is the newly added edge
coupling_cell = cell(p, q);

get_c(p, q);

dFrechet = CA(p, q);
coupling = get_coupling(p, q);

function coupling = get_coupling(i, j)
coupling = [];
uv_this = [i; j];
while(1)
tmp = coupling_cell{uv_this(1), uv_this(2)};
coupling = [tmp(:, 2) coupling];
uv_last = tmp(:, 1);
if (uv_last(1) == -1 )
break;
else
uv_this = uv_last;
end
end

end

function get_c(i, j)
if (CA(i, j) > -1) % value already available
return;
end

d = dist_fcn(P(:, i), Q(:, j));

if (i == 1 && j == 1)
CA(i, j) = d;
coupling_cell{i, j} = [-1 i; -1 j];
elseif (i > 1 && j == 1)
get_c(i - 1, j);
CA(i, j) = max([CA(i - 1, j), d]);
coupling_cell{i, j} = [i - 1, i; j, j];
elseif (i == 1 && j > 1)
get_c(i, j - 1);
CA(i, j) = max([CA(i, j - 1), d]);
coupling_cell{i, j} = [i, i; j - 1, j];
elseif (i > 1 && j > 1)
get_c(i - 1, j - 1);
get_c(i - 1, j);
get_c(i, j - 1);
[val, ind] = min([CA(i - 1, j - 1), CA(i - 1, j), CA(i, j - 1)]);
switch ind
case 1
coupling_cell{i, j} = [i - 1, i; j - 1, j];
case 2
coupling_cell{i, j} = [i - 1, i; j, j];
case 3
coupling_cell{i, j} = [i, i; j - 1, j];
end
CA(i, j) = max(val, d);
end

end

end