There are probably a bunch of ways to do this. Considering that the endpoint trajectories are not aligned to the reference grid, I figured I'd just do it this way, even if the error would have been negligible otherwise.
S = load('points.mat');
px0 = squeeze(S.px(:,:,1));
py0 = squeeze(S.py(:,:,1));
pxf = squeeze(S.px(:,:,end));
pyf = squeeze(S.py(:,:,end));
% get rid of nans
px0 = px0(~isnan(px0));
py0 = py0(~isnan(py0));
pxf = pxf(~isnan(pxf));
pyf = pyf(~isnan(pyf));
% order points to form a polygon
x = [px0; flipud(pxf)];
y = [py0; flipud(pyf)];
sh = polyshape(x,y);
plot(sh); grid on
axis tight
area(sh) % this is the area of the calculated polygon
Of course, the first dataset isn't as described by your image. If there's an offset or if px0, py0 are something other than the first vectors in in the arrays, then you'll have to make that change.
If you simply want to exclude the overlap area, then you can do this
% truncate overlap
x = [px0; max(flipud(pxf),interp1(py0,px0,pyf))];
sh = polyshape(x,y);
plot(sh); grid on
axis tight
area(sh)
If the overlapped regions are to be considered as some sort of negative area, then you can use a similar technique to calculate it (or just take the difference of the two methods given).
polyshape() will spit out warnings in either case, but they're just warnings. It's just telling you that the polygon intersects itself, which is easy enough to tell just by looking.
