That code was written so long ago, I forgot I ever wrote it. But also, it was written in the days when the convex hull and triangulation tools in MATLAB were far less mature/sophisticated than they are now.
If I look at your dataset, it is just a huge number of points that all fall on a nice integer lattice. You probably extracted them as pixels from an image.
plot(xy(:,1),xy(:,2),'.')
All of the codes in that toolbox generally first took the convex hull of the data. Anything inside the convex hull is meaningfless in terms of a bounding polygon anyway. It dramatically reduces the problem, since it needs only to work with the edges of the convex hull.
T = convhull(xy(:,1),xy(:,2));
plot(xy(T,1),xy(T,2),'-o')
The problem arises since your data lives purely on an integer lattice. And that means that at least a few of those edges were collinear edges in the simple convex hull. In turn, that got the code confused.
T2 = convhull(x,y,'Simplify',true);
plot(xy(T2,1),xy(T2,2),'r-s')
As you can see here, the simplify option was smart to replace those multiple collinear edges with a single edge. That resolves the problem in the code, and it will make the code run faster too.
I must post new versions of the codes in that toolbox.