Frechet Distance (discrete)

A scalar similarity measure between two curves. Points can be n-dimensional and custom distance metric is permitted.

Calculates the discrete Frechet distance between curve1 and curve2.

Simple usage:

frechetDist = frechetDistance(curve1, curve2)

Returned in frechetDist is the distcrete Frechet distance between the curves supplied in curve1 and curve2. The curves should have one point per row, and may be of arbitrary (but matching) dimensionality. In simple usage, the distance is the n-dimensional Euclidean distance, and endpoint flexibility is not allowed (i.e. the sequence of point couplings mustinclude pairing the starting points and pairing the ending points of each curve).

Complex usage:

[frechetDist, couplingSequence, frechetMatrix, distanceMatrix] = ...
   frechetDistance(curve1, curve2, distanceMatrixFunction, endpointFlexibilityFlag)

The output frechetDist and the inputs curve1 and curve2 are the same as in simple usage.

Optional input distanceMatrixFunction should be a function handle for a function which takes in curve1 and curve2 as inputs and which returns the pairwise distance matrix between the points of curve1 and the points of curve2. That is, if dm = distanceMatrixFunction(curve1, curve2), then dm(i,j) is the distance between the point curve1(i,:) and the point curve2(j,:). This allows the user to specify any arbitrary desired distance metric calculation, and would even allow a distance to be specified between points of differing dimensionalities (if such a thing were desired). If this input is empty ([]) or the number of arguments is less than 3, distanceMatrixFunction defaults to n-dimensional Euclidean distance (i.e. the L2 norm between points).

Optional input endpointFlexibilityFlag should be either true or false. If false, empty, or there are fewer than 4 input arguments, endpoint flexibility is not allowed, meaning that the sequence of point couplings must include pairing the starting points and ending points of each curve (which is standard for typical Frechet distance calculations). However, for some applications, this behavior is misleading about the dissimilarity between curves. For example, consider two periodic curves defined as follows:

x1 = linspace(0,2*pi,500)';
y1 = sin(10.*x1);
curve1 = [x1, y1];
x2 = x1;
y2 = sin(10.*(x2-0.1)); % shift 0.1 in x
curve2 = [x2, y2];
frechetDist = frechetDistance(curve1, curve2)
% = 0.8415
endpointFlexibilityFlag = true;
frechetDist = frechetDistance(curve1, curve2, [], endpointFlexibilityFlag)
% = 0.101

Intuitively, the Frechet distance between these curves should be 0.1, since curve2 is just curve1 shifted to the right by 0.1. However, the Frechet distance between these two curves comes out to 0.8415 rather than 0.1, because the standard processing forces the starting points to be paired together, and these points are 0.8415 apart. If the curves continued infinitely, then every point could be paired with one 0.1 away, but since the curves are out of phase and we are looking at the same range of x values, the endpoints are typically more than 0.1 apart. The endpointFlexibilityFlag allows this type of issue to be mitigated by relaxing the requirement that both starting points (and ending points) must be paired together. Instead, the distance between the first point of curve1 and the closest point on curve2 to it is found, and the distance between the first point of curve2 and the closest point on curve1 is found, and whichever of these is shorter forms the first link in the coupling sequence. Essentially, this allows ignoring the early part of one (but not both) of the two curves when calculating the Frechet distance. The same process is used at the other end of the curves, allowing the tail end of one of the curves to be ignored if there is a closer link between one of the ends and an earlier point on the other curve. If we run this function again on the curves defined above and allow endpoint flexibility, we find that the calculated Frechet distance falls to essentially 0.1 as expected (it's slightly greater because of the discrete sampling of the curves). It is recommended to think carefully about your data and whether it would be appropriate to allow endpoint flexibility or not. The easiest way to see the impact is to use plotCouplingSeq.m with a coupling sequence generated without endpoint flexibility and then with a coupling sequence generated with endpoint flexibility and compare them.

Output couplingSequence is an Nx2 matrix of point indices representing one possible sequence of curve1 to curve2 point couplings with the calculated Frechet distence; that is, the distance between the indicated pairs of points is equal to the Frechet distance for at least one pairing, and less than or equal to the Frechet distance for all other pairings. The first column of couplingSequence is a row index into curve1, and the second column of couplingSequence is a row index into curve2. Each column vector has the property that the next entry is either equal to or one greater than the previous entry, guaranteeing that curve points are visited in sequential order. If endpointFlexibility is not allowed, the first row of couplingSequence will always be [1,1] and the last row will be [n1,n2] if there are n1 points in curve1 and n2 points in curve2. If endpointFlexibility is allowed, then the first row will include at least one 1 and the last row will include at least one of n1 or n2. Note that there may be many possible coupling sequences with the same Frechet distance, since there may be many options for couplings which are not minimal but do not exceed the Frechet distance. This function tries to build a "good" coupling sequence by preferring couplings which are shorter distance to those which are longer, but even these are not guaranteed to be unique. The returned couplingSequence should be considered an illustrative example of a possible coupling which satisfies the calculated frechetDist.

Output frechetMatrix is, if endpointFlexibility is not allowed, an n1 by n2 matrix where the frechetMatrix(i,j) is the Frechet distance between the curves curve1(1:i,:) and curve2(1:j,:). Thus, the Frechet distance between the full curve1 and curve2 is frechetMatrix(n1,n2). The frechetMatrix is iteratively or recursively constructed starting from the pairwise distance matrix by considering the minimal longest link to get to a certain pairing of points. If endpointFlexibility is allowed, then frechetMatrix may be smaller than n1 by n2, because ignored points near either end will not be included. For example, if the first 3 points of curve1 are ignored and the last 2 points of curve2 are ignored, then the output frechetMatrix will be (n1-3) by (n2-2). The calculated frechetDist will still be the last-row, last-column entry of the frechetMatrix.

Output distanceMatrix is the full pairwise distance matrix between the points of curve1 and curve2. It will always be n1 by n2 regardless of whether endpoint flexibility is allowed.

Also included is plotCouplingSeq.m, which allows easy visualization of the results of frechetDistance processing. Links in the coupling sequence are shown as thin gray lines and any links which are at the Frechet distance are highlighted in green. This allows easy visual understanding of where the curves are most distant and also how the point pairings might be constrained. This function could also be easily modified to include things like customized link color or weight based on the pairing distance if desired.

This contribution was inspired by "Discrete Frechet Dist" and "Frechet Distance Calculator" entries on the FileExchange, after I discovered a bug in the Frechet Distance Calculator code (points could be visited in a non-monotonic order) and had trouble clearly understanding the code in Discrete Frechet Dist. My version includes multiple methods of calculation for the Frechet distance and is more extensively commented than either of the inspirations. My construction of the coupling sequence is inspired by, but an improvement on, the approach used in Discrete Frechet Dist. After noticing that coupling sequences produced by that approach often seemed awkward in the region after the Frechet link (because underconstrained), I added the heuristics to prefer shorter coupling links (an elastic leash in the leashed dog metaphor for Frechet processing) and, if link length is equal, to prefer more direct walks through the Frechet matrix. Finally, the ability to allow endpoint flexibility is entirely new and was inspired by the particular data I was evaluating (which was approximately periodic time series data).

Code to reproduce image:

%% To reproduce example image
t = linspace(0,4*pi,100);
x1 = t'; 
y1 = sin(x1);
x2 = t';
y2 = 1.5 .* sin(x2 + 0.3);

curve1 = [x1,y1];
curve2 = [x2,y2];

[fd,cs,fm,pd] = frechetDistance(curve1,curve2);
curveLabels = {'sin(x)', '1.5*sin(x+0.3)'};
plotCouplingSeq(curve1,curve2, cs, pd, curveLabels)
title(sprintf('Frechet distance = %0.2f',fd))