Alphashape boundary detection not working as desired

Given a complex set of (x, y) data I would like to create a hull that is not convex but not fully conforming. I (think I) prefer alphashape over boundary as the former lets me select a specific number of pixels for the radius. My final goal is to create a mask from the outline of the resulting shape. However, my attempt has failed.

Any; suggestions on how to make this work as intended?

load XY.mat

shp = alphaShape(x, y, 100, 'HoleThreshold', 99999999999999);
[~, BXY] = boundaryFacets(shp);

figure
scatter(x, y, '.k')
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

figure
plot(shp)
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

% This looks good.  I want everything green/black in the final mask.

figure
plot(BXY(:,1), BXY(:,2))
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

% This seems to have the correct points, but  there are lines both 
% inside and outside of the desired shape and for it to work right it 
% needs to be just the outside boarder of the above

figure
BW = poly2mask(BXY(:,1), BXY(:,2), max(x) + 500, max(y) + 500);
imshow(BW);

% Consequently, the mask looks all wrong.
% I want the mask to to be the entire area enclosed by the alphashape, 
% with no holes, and just the area enclosed by the alphashape.

[James Akula]

Closer. But still need help...

load XY.mat

shp = alphaShape(x, y, 100, 'HoleThreshold', 99999999999999);
[~, BXY] = boundaryFacets(shp);

subplot(2, 2, 1)
scatter(x, y, '.k')
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

subplot(2, 2, 2)
plot(shp)
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])
% This looks good

subplot(2, 2, 3)
plot(BXY(:,1), BXY(:,2))
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])
% This looks good, except there are crisscrossing lines and for it to work
% right it needs to be just the outside boarder

subplot(2, 2, 4)
[R, C] = find(ones([max(y) + 500, max(x)+500]));
tf = inShape(shp, C, R);
BW = false([max(y) + 500, max(x) + 500]);
BW(tf) = true;
imshow(BW)

% I was ablt to solve part of my problem--that is, the creation of a
% working mask.  But I would still like the border plot (lower left) to
% work, too!

Answer 1

[James Akula]

An answer that seems to work is to convert the alphashape into a triangulation and then use "boundaryshape."

load XY.mat

shp = alphaShape(x, y, 100, 'HoleThreshold', 99999999999999);

subplot(2, 2, 1)
scatter(x, y, '.k')
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

subplot(2, 2, 2)
plot(shp)
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

subplot(2, 2, 3)
[tri, P] = alphaTriangulation(shp);
TR = triangulation(tri, P);
polyshp = boundaryshape(TR);
plot(polyshp, 'EdgeColor', 'b', 'FaceColor', 'none')
axis ij equal off
axis([1, max(x) + 500, 1, max(y) + 500])

subplot(2, 2, 4)
ps_regions = regions(polyshp);
BW = false([max(y) + 500, max(x) + 500]);
for ii = 1:numel(ps_regions)
    BW = BW | poly2mask(...
        ps_regions(ii).Vertices(:,1), ps_regions(ii).Vertices(:,2), ...
        max(y) + 500, max(x) + 500);
end
imshow(BW)

Answer 2

[Image Analyst]

Try imfill

BW = imfill(BW, 'holes');

[James Akula]

Thank you for this advice, but filling the holes isn't the solution, as that would do "too much." The "outline" of the blue plot does not match that of the alphashape, as lines crisscross both internally and externally.

Answer 3

[Matt J]

Try bwlalphaclose() from this FEX submission,

https://www.mathworks.com/matlabcentral/fileexchange/90762-further-tools-for-analyzing-objects-in-n-dimensional-images

load XY
s=5;
x=round(x/s)+10; y=round(y/s); %cut down to more manageable sizes

shp = alphaShape(x,y, 100/s, 'HoleThreshold', 99999999999999/s);

Warning: Duplicate data points have been detected and removed.

BW = accumarray([y,x], true, [max(y) + 500/s, max(x) + 500/s]);

BW=bwlalphaclose(BW,shp);

imshow(BW,[]); shg
hold on
scatter(x,y,'r.')
hold off

[James Akula]

I did download this package but I struggled to make it work...I'm not sure that it actually does.

[Matt J]

Did you run the example in Examples.mlx?

[Matt J]

Here's what I get:

load XY
s=5;
x=round(x/s)+10; y=round(y/s); %cut down to more manageable sizes

shp = alphaShape(x,y, 100/s, 'HoleThreshold', 99999999999999/s);

Warning: Duplicate data points have been detected and removed.

BW = accumarray([y,x], true, [max(y) + 500/s, max(x) + 500/s]);

BW=bwlalphaclose(BW,shp);

imshow(BW,[]); shg
hold on
scatter(x,y,'r.')
hold off